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Abstract 

The  National  Air  Intelligence  Center,  WPAFB,  OH,  needs  to  predict  radar  beam 
refraction  with  greater  accuracy.  Hitherto,  beam  bending  has  been  predicted  using  four- 
thirds  earth  or  standard  atmosphere.  A  new  and  more  accurate  model  was  developed  for 
this  thesis  that  replaces  the  old  rules-of-thumb  with  a  mix  of  raytracing  and  climatology. 

Usually  a  microwave  radio  beam  traveling  through  the  atmosphere  bends  towards 
the  earth  with  a  radius  of  curvature  greater  than  the  earth’s  surface.  However,  seasonal 
and  climatic  variations  influence  the  amount  and  direction  of  bending,  and  at  times  create 
temperature  or  moisture  inversions  that  act  to  redirect  the  energy  along  the  earth’s  surface 
leaving  gaping  radio  holes  where  there  is  no  coverage. 

This  model  uses  iterative  raytracing  to  determine  the  most  direct  path  from  the 
radar  to  the  target  through  the  climatologically  predicted  refractive  atmosphere.  The 
amount  of  height  measurement  error  is  calculated  by  comparing  the  geographic  path  to 
the  refracted  path.  Only  vertical  refractivity  variation  is  taken  into  account,  and  the 
effects  of  ducting  and  exponential  refractivity  are  both  modeled. 

As  a  test,  the  model  computed  height  error  at  17  locations  worldwide  for  a  target 
at  10,000  feet  and  60  nautical  miles.  The  predicted  errors  varied  from  approximately  100 
feet  to  2260  feet  -  widely  varying  from  the  standard  atmosphere  predicted  height  error  of 
804  ft.  The  model  traces  to  all  targets  when  no  ducting  is  modeled,  to  all  targets  outside 
the  duct  with  surface  ducting,  and  to  some  targets  outside  the  duct  with  elevated  ducting, 
since  in  this  case  adjacent  rays  sometimes  cross,  causing  ambiguity  in  the  estimation. 
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A  CLIMATOLOGY-BASED  MODEL  FOR  STRATEGIC  PREDICTION  OF  RADAR 


BEAM  REFRACTION 


/.  Introduction 

As  early  as  1919  radio  scientists  were  investigating  the  effects  of  tropospheric  refraction 
or  bending,  of  radio  waves  (Kerr,  1951:2).  With  the  invention  of  radar,  accurate 
prediction  of  radio  beam  refraction  became  especially  important.  Since  one  assumes  a 
line-of-sight  path  to  determine  the  range  of  a  given  target  and  any  deviation  from  a 
straight  line  path  causes  erroneous  measurements.  Height  measurement  errors  are  of 
particular  concern  since  the  atmosphere  primarily  acts  to  bend  the  beam  in  a  vertical 
direction.  Beam  curvature  can  also  produce  range  measurement  errors,  but  they  are  less 
significant. 

With  the  proliferation  of  inexpensive  yet  powerful  computers,  radar  engineers 
turned  the  algorithms  compiled  by  people  like  Donald  Kerr,  Lamont  Blake,  Bean  and 
Dutton,  and  others  into  convenient  software  models  to  be  used  by  radar  engineers, 
evaluators,  and  mission  planners.  Recently,  the  National  Air  Intelligence  Center  (NAIC), 
at  Wright-Patterson  Air  Force  Base,  began  developing  AMBER,  a  general-purpose  radar 
range  prediction  model  they  intend  to  be  the  most  accurate  yet. 
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Tropospheric  Refraction 

Everyone  has  observed  the  apparent  bending  of  spoon  sitting  in  a  glass  of  water, 
or  the  shimmering  effect  of  the  air  over  tarmac  on  a  blisteringly  hot  day.  What  is  seen  in 
both  cases  results  from  the  bending  of  light  waves  as  they  pass  from  one  medium  to  a 
another  medium  with  different  indices  of  refraction.  This  phenomenon  not  only  applies 
to  light,  but  all  electromagnetic  wave  energy  including  the  ultra  high  and  higher  radio 
frequencies  used  in  modern  radar. 

The  index  of  refraction  of  air  depends  upon  temperature,  pressure,  and  moisture 
content  (humidity).  Thus,  when  a  radio  wave  propagates  through  the  troposphere,  a  layer 
of  atmosphere  which  extends  to  a  height  of  approximately  15  km,  it  experiences  a 
continually  changing  medium.  As  might  be  expected,  the  index  of  refraction  of  air  varies 
somewhat  predictably  with  altitude.  Furthermore,  it  is  also  governed  by  weather, 
geography,  time  of  day,  and  local  climate.  Tropospheric  refractivity  is  a  complex  and 
largely  unpredictable  quantity. 

For  better  or  for  worse,  it  is  exactly  that  quantity  which  the  radar  engineer  must 
predict  to  accurately  estimate  how  his  radar  beam  will  perform  at  a  given  time  and 
location.  Fortunately,  scientists  have  been  studying  the  troposphere  for  decades  and  have 
not  only  been  able  to  document  the  average  refractive  characteristics  of  the  atmosphere, 
but  have  done  extensive  analysis  on  the  climatological  variations  as  well.  For  example. 
Bean  and  Dutton’s  climatological  research  was  based  on  a  five-year  study  (Bean  and 
Dutton,  1966:109),  and  the  NCEP/NCAR  Reanalysis  Project  is  based  on  40  years  of 
measured  data  (Kalnay  et.  al,  1996). 
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As  these  and  other  studies  have  shown  over  and  over  again,  the  index  of  refraction 
on  an  average  decreases  exponentially  with  altitude.  The  gradient  of  this  decrease  again 
depends  upon  the  temperature,  pressure,  and  moisture  content  of  the  air.  Fortunately,  the 
air  at  the  surface  is  a  good  index  of  what  is  happening  in  the  upper  air  and  we  can 
construct  a  reasonable  approximation  to  the  actual  profile  using  the  surface  index  of 
refraction  and  the  initial  refractivity  gradient.  A  radar  beam  propagating  through  such  an 
atmosphere  will  gradually  bend  towards  the  earth  with  a  radius  of  curvature  much  greater 
than  that  of  the  earth  itself.  At  higher  initial  elevation  angles,  the  beam  bending  may  be 
insignificant,  but  at  angles  near  one  or  two  degrees  the  bending  can  be  extreme  causing 
height  measurement  errors  of  thousands  of  feet  at 'ranges  of  a  two  or  three  hundred 
nautical  miles.  Fortunately,  these  errors  can  be  predicted  with  some  accuracy. 

More  generally,  various  meteorological  and  climatic  effects  such  as  extreme  or 
rapid  heating  and  cooling,  sea  breezes,  thunderstorms  and  subsidence  will  cause  moisture 
and/or  temperature  inversions  creating  regions  of  superrefraction  in  which  the  radio 
waves  are  bent  so  much  they  are  funneled  along  the  earth’s  surface  for  large  distances. 
This  funneling  is  called  ducting,  and  can  occur  near  the  surface  as  well  as  at  higher 
altitudes  when  it  is  referred  to  as  elevated  ducting.  These  ducts  create  radio  holes, 
regions  above  (or  below)  the  duct  where  the  radio  wave  should  have  reached,  but  could 
not  because  it  was  redirected  away  from  its  expected  path.  All  these  effects  are  known  as 
anomalous  propagation,  though  in  many  places  a  troposphere  riddled  with  ducts  is  more 
common  than  the  smooth  monotonic  atmosphere  we  commonly  call  standard. 

Modeling  even  the  standard  (non-anomalous)  atmosphere  relies  on  either  current 
atmospheric  data  at  the  radar  site,  or,  if  that’s  not  available,  at  least  a  knowledge  of  the 
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local  climate  and  the  average  conditions  that  prevail.  The  anomalous  effects  can  also  be 
approximated,  but  with  less  certainty  since  they  change  continually. 

Problem  Statement 

The  current  model  used  by  NAIC  implements  the  basic  equations  found  in  Blake 
(1986: 179-188)  and  gives  the  option  to  use  either  an  exponential  refractivity  profile  based 
on  the  average  measurements  for  the  continental  United  States  or  a  low-altitude  linear 
approximation  to  that  known  as  the  effective  earth  approximation  or  4/3  earth.  Although 
these  approximations  are  better  than  no  refraction  model,  they  rarely  predict  the  true 
nature  of  the  atmosphere  with  any  accuracy. 

There  are  several  inherent  requirements  for  the  improved  model  desired  by  NAIC. 
First,  it  must  be  able  to  predict  the  radar’s  height  and  range  measurement  error 
performance  with  respect  to  a  given  target.  Second,  it  must  be  able  to  be  incorporated 
into  the  main  AMBER  model,  which  in  turn  requires  well-organized,  modular  code  with 
thorough  documentation,  and  specifically  requested,  prototyping  done  in  MATLAB®. 
Finally,  the  model  is  to  be  geared  towards  long  term  prediction.  That  is,  it  will  be  used 
for  hypothetical  scenarios,  or  scenarios  that  might  occur  some  time  in  the  future  without 
the  benefit  of  real-time  weather  data  from  the  radar  site.  However,  it  is  perfectly 
reasonable  to  expect  to  know  the  location  of  the  site,  the  time  of  year,  and  the  time  of  day 
the  radar  will  be  operating.  All  these  requirements  were  either  explicitly  or  implicitly 
specified  by  NAIC. 

This  third  requirement  suggests  a  model  based  on  climatology.  That  is,  the 
refractivity  profile  used  must  reflect  the  average  conditions  for  the  geographic  location  of 
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the  radar,  the  time  of  year,  and  the  time  of  day.  No  long-term  prediction  model  could 
make  use  of  day-to-day  (synoptic)  weather  phenomena  since  those  are  in  no  way 
predictable.  The  model  should,  if  possible,  be  able  to  take  into  account  both  the  smooth, 
monotonic  characteristics  of  the  standard  atmosphere  and  the  anomalous  propagation 
conditions  associated  with  ducting,  assuming  enough  climatological  knowledge  is 
available.  At  first,  one  might  dismiss  the  possibility  of  being  able  to  predict  ducting 
using  climatological  information,  but  we  can  at  least  draw  statistics  from  historical  data 
with  respect  to  the  chances  of  a  duct  occurring  and  its  average  characteristics.  So  ducting 
also  should  be  taken  into  account. 

Of  course,  this  is  not  the  first  time  atmospheric  refraction  has  been  modeled  using 
data  other  than  the  standard  atmosphere.  The  Naval  Command  Control  and  Ocean 
Surveillance  Center  (NCCOSC)  in  San  Diego  in  1976  (Hitney  and  Richter,  1976) 
described  the  first  of  their  models,  Integrated  Refractive  Effects  Prediction  System 
(IREPS).  The  upgrades,  Engineering  Refractive  Effects  Prediction  System  (EREPS) 
(Patterson,  1994)  and,  most  recently.  Advanced  Refractive  Effects  Prediction  System 
(AREPS)  (Patterson,  1998)  have  followed.  AREPS  uses  a  combination  of  geometrical 
optics  (ray-tracing)  and  a  model  known  as  the  electromagnetic  parabolic  equation  to 
calculate  power  loss  across  a  range-height  field  (Hitney  &  Richter,  1976,  Patterson, 
1998).  AREPS  accepts  a  variety  of  data  input  methods  including  direct  user  data  entry, 
and  an  extensive  (370  stations  worldwide)  climatological  database  put  together  by  the 
Navy  in  1987  (Patterson,  1987).  In  1982,  Abel,  et.  al,  described  their  model  RAYTRA 
(short  for  Raytrace),  a  comprehensive  raytracing  model  which  is  particularly  useful  for 
point-to-point  calculation,  unlike  the  others  which  provide  loss  predictions  over  a  whole 
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range  of  ranges  and  heights.  Another  parabolic  equation  model,  the  Tropospheric 
Electromagnetic  Parabolic  Equation  Routine  (TEMPER),  developed  by  Johns  Hopkins 
University  Applied  Physics,  and  based  on  the  work  of  Ko,  Sari  and  Skura  (1983),  some 
advanced  techniques  to  more  accurately  model  anomalous  conditions  in  the  atmosphere. 
Another  model  developed  at  the  NCCOSC  (then  the  Naval  Ocean  Systems  Center)  is 
VTRPE,  a  third  parabolic  wave  equation-based  model  developed  by  Frank  Ryan  (1991). 
All  of  these  models  have  been  tried  and  tested  and  each,  to  a  more  or  less  extent,  is 
available  for  use.  Of  course,  none  is  designed  to  be  integrated  into  AMBER  and  only 
RAYTRA  is  wholly  concerned  with  point-to-point  height  and  range  calculations. 

The  problem,  then,  is  to  develop  a  customized  model  of  atmospheric  refraction 
which  will  calculate  the  height  and  range  measurement  errors  for  a  given  radar  location, 
time  of  year,  time  of  day,  elevation,  and  a  given  target  height,  range,  and  azimuthal 
bearing  from  the  radar.  The  model  should  rely  on  climatologically-based  refractivity 
data,  and  use  either  an  existing  model/algorithm  or,  if  necessary,  develop  a  new  one. 
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n.  Literature  Review 


The  nature  of  microwave  refractivity  and  its  governing  atmospherics  is  a  complex 
subject.  As  with  any  real  system,  the  factors  affecting  the  turbulent  troposphere  must  be 
prioritized  and  studied  according  to  the  amount  of  influence  each  has,  using  idealizations 
where  necessary.  Hence,  this  chapter  can  not  be  exhaustive.  Instead,  it  highlights 
information  appropriate  to  the  modeling  presented  in  the  next  chapter.  First,  simple 
refractivity  and  anomalous  refractive  behavior  is  described.  Next,  the  effect  of 
refractivity  on  propagation  of  microwave  radiation  is  presented,  including  some  of  the 
raytracing  mathematics  used  to  predict  beam  bending.  Following  that  is  a  discussion  of 
the  major  meteorological  and  climatic  factors  that  influence  refractivity.  Finally,  to 
illustrate  the  effect  of  these  factors,  a  brief  survey  of  world  climates  is  included. 

N-Refractivity 

A  canoe  paddle  is  dipped  into  the  water  and  it  seems  to  bend.  The  thirsty  desert 
traveler  peers  vainly  at  the  pools  of  water  on  the  horizon,  only  to  find  he  has  been  looking 
at  a  mirage.  The  distortion  of  light  that  is  observed  in  both  cases  is  caused  by  a  gradient, 
either  discrete  (as  in  the  former  case  above)  or  continuous  (as  in  the  latter),  in  the 
electromagnetic  properties  of  the  media  through  which  the  light  travels.  These 
electromagnetic  properties,  specifically,  the  dielectric  constant,  8,  and  permeability,  |x, 
affect  the  speed  at  which  electromagnetic  (E-M)  waves  propagate.  To  wit,  the  ratio  of 
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the  speed  of  the  electromagnetic  wave  in  a  vacuum  to  the  speed  of  electromagnetic 
energy  in  the  medium  is  given  by 


V 


where  n  is  commonly  known  as  the  index  of  refraction.  Since  n  in  air  is  approximately 
1.0003,  a  more  convenient  representation  has  been  adopted,  that  is 

N  =  in-\)*l(f  ,  (2-2) 

where  N  is  called  the  refractivity  of  the  medium. 

The  simplest  way  to  understand  how  a  refractivity  gradient  causes  an  E-M  wave  to 
bend  is  by  considering  an  interface  between  two  media  of  differing  refractivity.  As  the 
wavefront  approaches  the  interface  at  some  oblique  angle,  the  portion  of  the  wavefront 
that  reaches  the  second  medium  first  will  either  speed  up  or  slow  down  depending  on 
which  medium  has  greater  refractivity.  This  speed  differential  in  the  wavefront  will  cause 
it  to  pivot  at  the  interface,  much  like  a  tractor  does  when  the  brake  is  applied  on  one  side 
or  the  other.  The  result  is  an  apparent  bending  of  the  wave  as  it  enters  the  new  medium. 
This  bending  is  governed  by  Snell’s  Law  of  Refraction  which  is  put  to  direct  use  in  the 
later  section  on  raytracing. 

The  earth’s  atmosphere,  of  course,  does  not  consist  of  two  homogeneous  media, 
but  of  a  mixture  of  various  gases  with  combined  pressure,  temperature,  and  water  vapor 
content  that  change  with  time,  altitude  and  horizontal  distance.  To  study  radio  beam 
bending  in  some  local  area,  however,  it  is  reasonable  to  ignore  the  horizontal  variation  of 
these  atmospheric  properties  and  consider  only  how  they  vary  with  height  at  a  given  time. 
This  assumption  is  convenient  and,  in  truth,  necessary  to  limit  the  complexity  of  the 
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problem  to  within  practical  bounds  (Kerr,  1951:46).  As  it  turns  out,  the  refractivity  of  air 
is  defined  by  (Kerr,1951:  13;  Bean  and  Dutton,  1966:  7) 
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where  T  is  temperature  in  Kelvins,  p  is  pressure  in  millibars,  and  e  is  the  partial  pressure 
of  water  vapor  in  millibars.  Although  sources  differ  somewhat  with  regard  to  the  value  of 
the  constants,  these,  given  by  Smith  and  Weintraub,  are  in  widespread  use  (Bean  and 
Dutton,  1966:  7,  Smith  and  Weintraub,  1953).  Further,  this  formulation  ofiVis  accurate 
for  frequencies  up  to  at  least  24  gigahertz  (Bean,  6).  Since  most,  if  not  all,  ground-based 
radars  operate  between  the  HF  (3-30  megahertz)  and  C  (4-8  gigahertz)  bands,  and  most 
airborne  radars  in  the  X  (8-12.5  gigahertz)  and  Ku  (12.5-18  gigahertz)  bands  (84th  Radar 
Evaluation  Squadron,  1998:  83-87),  this  limitation  is  acceptable.  Caution  may  be 
necessary  when  working  with  some  of  the  newest  ground  search  and  terrain  avoidance 
applications  for  which  Ka  (26.5-40  gigahertz)  and  millimeter-wave  bands  are  being 
exploited. 

Using  a  small  weather  set  called  a  radiosonde  that  is  carried  aloft  by  balloon, 
scientists  measure  pressure,  temperature,  and  water  vapor  at  any  desired  levels  up  through 
the  atmosphere  to  precisely  determine  the  actual  profile.  Radiosonde  measurements  have 
shown  refractivity  generally  decreases  exponentially  with  altitude.  Such  a  gradient  causes 
a  radio  wave  to  bend  downwards,  but  with  a  radius  of  curvature  much  greater  than  that  of 
the  earth’s  surface,  eliminating  the  possibility  of  the  wave  returning  to  the  earth.  Typical 
behavior  of  the  refractivity  profile  can  be  observed  readily  in  the  sample  data  for  Buffalo, 


NY  shown  in  Figure  2. 1 . 
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Figure  2.1  Sample  iV  Profile  Calculated  from  Radiosonde  Data  in  the  National  Climatic 
Data  Center  (NOAA)  Database. 


Several  techniques  have  been  developed  to  mathematically  approximate  this 
standard  refractivity  profile.  The  traditional  method  used  by  radar  engineers  for  quick 
calculations  is  the  effective  earth  radius,  or  4/3-earth  approximation.  This  model, 
originally  proposed  by  Schelleng,  Burrows,  and  Ferrell  (1933)  is  obtained  by 
straightening  the  curved  path  of  the  radio  waves  without  changing  the  height  of  the 
wavefront  at  any  point  along  the  path.  The  effect  of  this  distortion  is  to  also  partially 
flatten  the  earth’s  surface  until  it  has  a  larger  radius  than  actual  (see  Figure  2.2).  The 
effective  earth  radius  is  defined  by 

=  ka  (2-4) 

where  a  is  the  actual  earth  radius  and 
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where  0  is  the  takeoff  angle  of  the  radar  beam.  Setting 


reality  as  shown  in  Figure  2.3. 


Actual  radio  ray 
Effective  radio  ray 


Figure  2.2  Effective  Earth  Approximation 


Buffalo,  NY  --  1  Jan  92  --  a.m. 


Figure  2.3  Three  standard  models  overlaying  measured  data:  4/3  Effective  Earth  Radius 
model,  Single  Exponential  Model  with  gradient  determined  using  CRPL  Standard 
Atmosphere,  and  Single  Exponential  Model  with  gradient  drawn  from  climatological 
studies 


The  next  model  more  closely  approaches  the  exponential  nature  of  the  A-profile. 
In  this  one,  the  surface  refractivity.  As,  and  the  initial  gradient  up  to  1000  meters  above 
the  surface,  AN,  are  used  to  calculate  a  single  exponential  distribution  of  A: 


where 


N  =  Ns  &xp{-c^(h-hs)}. 


,  ,  Ns 

c,  =  In — ^  =  In- 


(2-8) 


A,^  A,+AA 


(2-9) 


Values  for  As  may  be  readily  obtained  from  surface  measurements  or  from  climatology 
studies,  that  is  studies  that  investigate  the  average  values  for  a  particular  location  during  a 
particular  time  of  year.  Lacking  this  sort  of  information,  one  can  simply  use  the  value 
generally  considered  to  represent  average  conditions  in  the  U.S.,  As=313  (This  would 
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correspond  to  k=4/3  in  the  linear  model).  Similarly,  AN  may  be  derived  from  radiosonde 
measurements,  or  in  some  cases,  climatology;  however,  this  value  may  also  be  calculated 
from  Ns  using. 


-AN  =  132Qxp{Q.005571Ns}  ■  (2-10) 

Equations  (2-8),  (2-9),  and  (2-10)  are  collectively  known  as  the  Central  Radio 
Propagation  Laboratory  Exponential  Reference  Atmosphere  (Bean  and  Thayer,  1959). 
Bean  and  Thayer  obtained  this  last  equation  after  analyzing  6-year  means  for  45  U.S. 
weather  stations  representing  all  kinds  of  climatic  conditions.  Though  not  as  accurate  as 
a  local,  seasonal  mean,  this  relationship  is  a  good  approximation  and  is  handy  when  no 
data  is  available  for  the  upper  air.  Figure  2.3  illustrates  the  exponential  profiles,  one 
using  CRPL  Reference  Atmosphere  and  the  other  using  climatology  averages  for  the  local 
area  and  time  of  year. 

Further  accuracy  in  modeling  standard  refractivity  may  be  accomplished  using  the 
bi-exponential,  or  as  it  is  commonly  known  today,  the  tri-exponential  model,  described 
thoroughly  in  Bean  and  Dutton  (1966:  311-322).  Refractivity,  N,  can  be  considered  to  be 
composed  of  a  dry  term. 


D  = 


11.6P 

T 


(2-11) 


and  a  wet  term. 


which  combine  to  give. 


3.73xl0^e 

t2 


(2-12) 
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N(h)  =  DgCxp 


(2-13) 


h  [  \  h 

—  f  +  WoCxpi  - 


H. 
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where  Do  and  Wo  are  the  dry  and  wet  terms,  respectively,  at  the  surface;  h  is  the  height 
above  the  surface;  and  Ha  and  Hw  are  the  scale  heights  of  D  and  W,  respectively.  The 
scale  height  is,  in  each  case,  the  height  at  which  the  value  of  the  atmospheric  property  has 
decreased  to  1/e  of  its  surface  value.  The  primary  advantage  of  the  tri-exponential  model 
is  that  the  contributions  of  the  partial  pressure  of  water  vapor  content  and  the  dry  pressure 
can  be  examined  separately  to  obtain  a  clearer  understanding  of  the  makeup  of  the 
atmosphere.  To  illustrate.  Bean  and  Dutton  (1966;  312)  included  the  following  table, 
which  includes  typical  values  for  three  distinct  climate  types: 


Table  2.  1  T5^ical  average  values  of  the  dry  and  wet  components  of  N 


Station  and  Climate 

Do 

Wo 

No 

Isachsen  t78°  50’  N).  arctic . 

Washington,  D.C.  (38°  50’  N), 

Canton  Island  (2°  46’  S),  tropic . 

332.0 

266.1 

259.4 

0.8 

58.5 

111.9 

332.8 

324.6 

371.3 

In  the  cold,  dry  arctic,  the  dry  component  makes  up  the  vast  majority  of  No,  but  in  the 
tropics,  where  humidity  is  high,  the  wet  term  makes  up  a  much  greater  portion  of  the  total 
refractivity.  Hence,  this  model  is  particularly  useful  for  climatological  studies  of  N. 

The  preceding  models  are  predicated  upon  the  atmosphere  having  a  smooth,  well- 
behaved,  monotonically  decreasing  index  of  refraction.  Though  it  may  be  tempting  to 
regard  this  sort  of  atmospheric  behavior  as  “normal,”  use  of  this  particular  term  has  been 
avoided  so  far  for  the  simple  reason  that  standard  refraction  is,  in  many  cases,  not  the 
norm.  There  are  many  places  and  seasons  during  which  refractive  anomalies  have  a  much 
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greater  influence  upon  the  bending  of  radio  waves  than  does  the  pure  exponential 
distribution. 

To  more  easily  observe  anomalous  behavior,  it  is  necessary  to  introduce  a  new 
term,  the  modified  index  of  refraction, 

M  =  N{h)  +  ih/aJlO^  (2-14) 

with  a  gradient, 

dM  dN  10^ 

=  + -  (2-15) 

dn  dh 

where  h  is  height  and  a®  is  the  radius  of  the  earth.  To  understand  the  physical  significance 
and  usefulness  of  this  parameter,  consider  that  the  M-gradient,  dM/dh,  goes  to  zero  at  any 
altitude  at  which  a  wave  launched  horizontally  travels  a  curved  path  that  is  concentric 
with  the  surface  of  the  earth.  In  a  standard  atmosphere,  M  will  increase  monotonically 
and  smoothly  as  shown  in  the  modified  version  (Figure  2.4)  of  the  now  familiar  Buffalo 
profile.  More  generally,  however,  meteorological  and  climatic  conditions  will  combine 
to  alter  this  simple,  pleasing  contour. 

Though  the  conditions  that  cause  refractive  anomalies  are  discussed  in  a  later 
section,  it  is  appropriate  here  to  understand  the  fundamental  t5'pes  of  irregularity.  Kerr 
(1951:  14,15)  classifies  non-standard  refractivity  characteristics  into  two  categories:  A 
layer  in  which  the  M-gradient  is  greater  than  standard  throughout  is  termed  substandard 
because  radio  frequency  propagation  in  a  sufficiently  thick  substandard  layer  is  usually 
poorer  than  expected.  Similarly,  a  layer  with  gradient  less  than  standard  is  called 
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superstandard,  because  propagation  performance  in  a  deep  enough  region  of  this  type  is 
generally  better  than  expected.  Note  that  whether  a  layer  has  an  M-gradient  between  zero 
and  standard,  or  an  M-gradient  less  than  0,  it  is  always  labeled  superstandard;  however, 
the  condition  leading  to  a  negative  M-gradient  is  more  strictly  termed  an  inversion. 

Figure  2.5  illustrates  idealized  M-profiles  combining  the  standard  profile  with  the  chief 
anomalies. 
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Figure  2.4  Sample  N  Profile  Calculated  from  Radiosonde  Data  in  the  National  Climatic 
Data  Center  (NO  A  A)  Database. 

Of  course,  this  description  of  N-  and  M-profiles  is  only  half  the  story.  Of  central 
importance  to  this  thesis  is  exactly  how  microwave  radiation  is  refracted  by  these  various 
conditions  in  the  troposphere. 
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Modified  Index,  M 


Figure  2.5  Idealized  modified  index  profiles:  (A)  Substandard  surface  layer;  (B)  profile 
for  standard  refraction;  (C)  superstandard  surface  layer;  (D)  superstandard  surface  layer 
with  surface  duct;  (E)  elevated  superstandard  surface  layer  with  surface  duct;  (F)  elevated 
superstandard  layer  with  elevated  ducts;  (G)  surface  and  elevated  superstandard  layers 
with  both  surface  and  elevated  ducts.  In  all  cases  the  duct  extends  from  aiob  and  from 
a’iob’.  (Kerr,  1951:  14) 


Tropospheric  Refraction  of  Microwave  Radio  Waves 

In  free  space  electromagnetic  energy  travels  in  a  straight  line.  In  fact,  the  same  is 
true  for  any  perfectly  homogeneous  medium.  However,  when  the  electrical  and  magnetic 
properties  along  the  path  begin  to  change,  segments  of  the  wavefront  begin  to  travel  at 
different  speeds  causing  the  wavefront  to  change  direction. 

Consider,  as  the  simplest  case,  a  wavefront  approaching  at  an  oblique  angle  the 
boundary  between  two  homogeneous  media  as  in  Figure  2.6.  Snell’s  law, 

n,  n, 

— ^  =  (2-16; 

cosa2  cosa, 

predicts  that  the  wavefront  will  turn  towards  the  normal  when  passing  into  a  medium 
with  higher  index  of  refraction.  This  is  what  happens  in  the  case  of  a  single,  discrete 
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change  in  n.  If  the  index  of  refraction  continues  to  change  stepwise  at  equally  spaced 
intervals,  the  wavefront  will  undergo  a  series  of  direction  changes.  Further,  if  the  spacing 
between  steps  is  reduced  to  a  differential  distance  and  the  index  of  refraction  step  turned 
into  a  gradient,  the  path  the  wave  follows  will  become  a  smooth  curve.  This  is  what 
happens  in  the  troposphere. 


Figure  2.6  Simplest  Case  of  Refraction 

The  different  classifications  of  refractivity  gradient  discussed  in  the  previous 
section  lead  to  different  wave  propagation  behaviors  as  shown  in  Figure  2.7.  A 
substandard  refractivity  gradient  causes  the  electromagnetic  wave  path  to  bend  upwards. 
Standard  refraction  causes  results  in  a  path  that  bends  down,  but  with  a  radius  of 
curvature  considerably  greater  than  that  of  the  earth.  Hence  the  beam  will  not  return  to 
earth.  Similarly,  superstandard  refraction  will  cause  the  beam  to  bend  down,  but  with  a 
smaller  radius  of  curvature.  These  three  types  of  bending  are  shown  again  in  Figure  2.8, 
except  that  here  the  earth  has  been  flattened,  changing  the  apparent  curve  of  the  rays. 
When  the  radius  of  curvature  is  small  enough  to  cause  the  elevation  angle  of  the 
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Figure  2.7  Four  categories  of  refractive  behavior:  (a)  subrefraction,  (b)  standard 
refraction,  (c)  superrefraction,  (d)  superrefraction  with  ducting 


Figure  2.8  (a)  1.  Substandard,  2.  Standard,  3.  Superstandard  profiles;  (b)  1.  Substandard, 
2.  Standard,  3.  Superstandard  bending 

beam  to  become  zero  or  negative  with  respect  to  the  local  tangent,  the  anomalous  form  of 
propagation  know  as  trapping,  or  ducting,  can  occur.  A  radio  beam  trapped  in  a  duct  can 
be  funneled  along  parallel  to  the  earth’s  surface  for  great  distances,  often  farther  than  it 
would  otherwise  be  able  to  travel  under  normal  line-of-sight  conditions.  This  said, 
ducting  can  occur  at  various  altitudes  and  for  a  variety  of  reasons.  Fundamentally,  a 
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particular  duct  can  be  classed  as  one  of  two  types:  If  it  occurs  because  of  an  inversion  at 
the  surface,  it  is  called  a  surface  duct,  or  ground-based  duct,  if  it  is  produced  by  an 
inversion  some  distance  above  the  earth,  it  is  called  an  elevated  duct  (see  Figure  2.5). 

Each  type  is  constructed  and  behaves  differently  from  the  other. 

Recall  from  the  previous  section  that  a  horizontal  ray  path  is  caused  by  an 
vanishing  M-gradient.  If  the  M-gradient  goes  negative,  the  ray  will  be  bent  down  towards 
the  earth.  When  the  initial  slope  (i.e.  near  the  surface)  of  the  M-profile  is  negative,  the 
result  is  a  surface  duct  (see  Figure  2.9).  If  the  initial  elevation  angle  of  a  particular  ray 
launched  from  ground  level  is  low  enough,  the  beam  will  be  bent  back  to  the  earth. 
Depending  on  the  electromagnetic  characteristics  of  the  ground,  the  beam  may  be 
reflected  only  to  be  bent  back  and  reflected  again.  In  this  way,  the  energy  can  be  funneled 
along  the  surface  for  great  distances.  Field  theory  can  account  for  this  behavior  using 
principles  similar  to  those  used  in  the  analysis  of  waveguides.  Kerr  (1951:  18-21) 
determined,  using  waveguide  theory,  the  maximum  wavelength  able  to  be  trapped  by  a 
surface  duct.  It  turns  out  the  limiting  factor  is  duct  height: 

X^^^=0mAd^’\  (2-17) 

where  'Kmax  is  the  maximum  trappable  wavelength  in  centimeters,  and  d,  in  feet,  is  the 
thickness  of  the  layer  with  negative  dM/dh.  Hence,  any  surface  duct  higher  than  500  feet 
will  trap  most  microwave  frequencies. 

If  the  negative  M-gradient  occurs  at  some  altitude  above  the  surface,  an  elevated 
duct  is  formed  (see  Figure  2.10).  For  the  most  part,  only  those  rays  originating  inside  the 
duct  may  be  trapped.  The  only  exception  is  when  a  ray  is  launched  above  the  top  of  the 
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duct  at  a  critical  angle  below  the  horizontal.  At  this  angle,  the  ray  can  be  trapped  by  a 
very  small,  hence  unstable,  trapping  layer  at  the  top  of  the  duct.  Any  other  ray  entering 
from  the  outside  will  only  have  its  course  altered  by  the  duct.  Rays  launched  from  inside 
the  duct  may  or  may  not  be  trapped,  again  depending  on  their  elevation  angle.  A  ray  with 
a  steep  enough  elevation  angle  will  not  be  bent  sufficiently  to  be  trapped  before  it  escapes 
the  duct  (Livingston,  1970:  105-114). 


(a)  (b)  (c) 


Figure  2.9  Surface  duct,  (a)  typical  profile,  (b)  paths  of  rays  launched  into  the  duct  from 
the  ground,  (c)  paths  of  rays  launched  from  above  the  top  of  the  duct  (adapted  from 
Livingston,  1970:  109) 

Although  the  calculations  involved  in  modeling  ducting  involve  special 
techniques,  we  may  calculate  the  tamer  refractive  behavior  with  considerable  accuracy 
using  the  techniques  of  geometric  optics,  more  commonly  known  as  raytracing  (Kerr, 
1951:  41).  A  sister  science,  physical  optics,  based  on  Maxwell’s  equations,  is  also  useful, 
but  is  more 
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Figure  2.10.  Elevated  duct,  (a)  profile,  (b)  paths  of  rays  launched  at  various  angles  both 
within  and  outside  the  duct  (adapted  from  Livingston,  1970:  111-114) 


often  employed  to  determine  the  power  losses  associated  with  refraction.  Since  the  goal 
of  this  thesis  is  to  more  accurately  predict  height  and  range  errors,  the  geometrical 
approach  will  suffice  as  long  as  certain  limits  are  kept  in  mind  and  taken  into  account. 
First,  the  refractive  index  must  not  vary  appreciably  in  one  wavelength;  and  second, 
neighboring  rays  must  remain  close  to  parallel  within  one  wavelength  (Kerr,  1951:  54), 
i.e.  they  cannot  cross.  It  follows  that  these  methods  become  more  and  more  accurate  as 
frequency  increases.  Additionally,  to  properly  use  Snell’s  Law,  we  must  consider  the 
atmosphere  as  a  series  of  differentially  thin  spherical  shells-a  stratified  model;  and  the  N- 
gradient  must  only  vary  with  height.  Assuming  these  limits  are  satisfied,  we  can  use  the 
following  form  of  Snell’s  Law  of  Refraction  (adapted  from  Kerr,  1951:  48-49): 

(a+li^)ni  costti  =(a  +  h2)n^  cosa^  (2-18) 
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where  a  is  the  earth’s  radius,  hj  and  /i2  are  the  upper  and  lower  boundaries  of  a  given 
shell,  H]  and  n2  are  the  indices  of  refraction  at  the  shell  boundaries,  and  ai  and  a2  are  the 
elevation  angles  of  the  ray  at  the  shell  boundaries  (Figure  2.1 1). 

To  calculate  the  bending,  and  more  specifically  (for  the  height/range  error 
application)  the  location  of  the  endpoint  of  the  ray,  we  must  sum  the  contributions  of  all 
the  shells  to  the  overall  bending.  So,  knowing  the  1)  initial  takeoff  angle,  2)  the  height 
increment,  and  3)  the  index  of  refraction  as  a  function  of  height,  we  begin  calculating  the 


Figure  2. 1 1.  Fundamental  raytracing  geometry  for  refraction  through  a  single  shell 
(adapted  from  Bean  and  Dutton,  1966:  50) 

effect  of  the  first  layer.  First,  we  use  Snell’s  Law  (2-17)  to  solve  for  a2,  the  elevation 
angle  at  the  upper  boundary  of  the  layer.  Next  we  calculate  the  total  bending,  y/i,  through 
the  layer,  using  a  relationship  derived  by  Weisbrod  and  Anderson  (1958;  Abel,  et.  ai, 
982:  16), 
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(2-19) 


2(n,  -n.)xlO^ 

^ - 

tana, -I- tan a2 

Finally,  the  subtense,  P,,  is  calculated  using  the  following  formula  (Abel,  et.  al,  1982: 

16), 

Pi  =\|/, +a2-ai,  (2-20) 

which  can  be  derived  using  the  quadrilateral  defined  by  the  angles  p,  ai  +  90°,  180°  -  y, 
and  90°  -  a2.  Then  we  simply  sum  the  subtenses, 

L 

Ktal  =  SP/  (2-21) 

l=\ 

where  L  is  the  total  number  of  layers,  to  find  the  subtense  between  the  ray  starting  and 
ending  points.  This  value,  along  with  the  height  of  the  ray  ending  point  and  the  Law  of 
Cosines,  can  be  used  to  calculate  the  geometric  (actual)  range,  Rg,  from  the  ray  starting 
point  to  its  ending  point  (Abel,  et.  al,  1982:  16): 

Rg  =  +  {a +  hff'  -2{a  +  hXa+  )  cos  P  ,  (2-22) 

where  hj  and  hf  are  the  initial  and  final  heights,  respectively.  Further,  the  length  of  the 
path  itself,  known  as  the  apparent  range,  Ra,  can  be  approximated  using  the  individual 
layer  subtenses,  and  boundary  heights  to  calculate  the  geometric  range  traversed  within  a 
layer,  and  summing: 

L  _ 

K  =  X'v(^‘‘‘^/)^  ■^(^  +  ^/+i)^  “2(^  +  ^/)(^  +  ^/+i)cosp,  (2-23) 

/=i 

These  calculations  may  be  applied  consistently  with  the  raytracing  limitations  defined 
above.  It  should  be  noted  that  an  alternate  and,  in  fact,  more  succinct  formula  for  Ra  is 
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derived  in  Blake  (1986:  180-182)  and  other  places.  The  Abel  method,  however,  better 
lends  itself  to  the  estimation  and  iteration  process  used  to  do  point-to-point  tracing. 

Using  these  techniques  in  the  presence  of  ducting,  the  most  interesting  refraction 
anomaly,  would  seem  not  to  be  sanctioned  because  a  refractivity  inversion  generally 
violates  the  first  raytracing  condition,  that  is,  the  gradient  may  not  change  appreciably  in 
one  wavelength.  However,  the  Navy,  in  the  RAYS  routine  of  their  EREPS  model 
(Patterson,  1994:  1 18-121),  applied  raytracing  to  the  duct  problem.  The  way  they 
constructed  the  ducting  profiles,  along  with  further  development  of  the  methods  used  to 
model  standard  refractivity,  are  described  and  applied  in  chapter  3. 

In  order  to  predict  the  refractivity  profiles  that  are  the  core  element  of  the  raytrace 
model,  it  is  necessary  to  understand  what  affects  tropospheric  refractivity  gradient  and 
how. 

Meteorological  Phenomena  that  Affect  the  Refractivity  Profile* 

To  understand  the  various  phenomena  that  control  tropospheric  refractivity  it  is 
necessary  to  establish  reference,  or  baseline  conditions  to  which  all  other  circumstances 
can  be  compared.  For  this  purpose,  scientists  have  defined  a  column  of  air  that  is 
completely  mixed,  or  stirred,  as  homogeneous.  In  such  vertically  homogeneous  air,  the 
vertical  variation  in  temperature  is  due  solely  to  the  change  in  pressure,  and  water  vapor 
concentration  is  independent  of  height.  Specifically,  the  decrease  in  temperature  with 
altitude,  or  lapse  rate,  is  adiabatic;  that  is,  it  represents  the  temperature  decrease  of  air 

*  The  material  in  this  section,  unless  otherwise  marked,  is  taken  from  Kerr,  1951:  181-293. 
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which  rises  and  cools  adiabatically.  Although  water  vapor  content  is  independent  of 
height,  the  vapor  pressure,  dew-point,  and  wet-bulb  temperature  do  vary  with  height 
because  they  depend  upon  changes  in  the  total  pressure  of  air.  Because  the  temperature 
lapse  rate  is  adiabatic,  the  air  will  not  move  vertically  unless  it  is  acted  upon  by  some 
external  influence.  Hence,  vertically  homogeneous  air  is  said  to  be  in  neutral 
equilibrium.  Likewise,  an  air  column  with  temperature  that  decreases  with  height  more 
slowly  than  the  adiabatic  lapse  rate  is  said  to  be  in  stable  equilibrium  and  an  air  column 
with  temperature  that  decreases  faster  than  adiabatic  is  said  to  be  in  unstable  equilibrium. 
The  vertically  homogeneous  air  standard  is  particularly  useful  because  it  is  the  only 
simple  distribution  that  occurs  often,  particularly  at  lower  altitudes.  When  it  does  occur, 
the  layer  is  bounded  on  the  top  by  the  altitude  at  which  condensation  occurs  (cloud 
height— on  the  order  of  10,000  feet),  and  on  the  bottom  by  the  top  of  the  surface  layer  of 
air  in  which  moisture  and  heat  are  exchanged  with  the  ground  (approximately  50  feet). 
Furthermore,  its  refractivity  characteristics  are  close  to  the  standard  atmosphere.  Drawing 
an  analogy  to  4/3  earth,  well-mixed  air  would  be  approximately  6/5  (or  3.6/3)  earth. 

To  arrive  at  this  homogeneous  state,  air  is  mixed  by  three  main  processes, 
convection,  eddy  turbulence,  and  molecular  diffusion.  Convection,  the  most  broad¬ 
ranging  process,  occurs  when  there  is  a  heat  source  increasing  the  temperature  of  the  air 
at  the  bottom  of  the  column.  When  this  happens,  the  heated  air  expands,  becoming  more 
buoyant,  and  begins  to  rise.  Cooler  air  then  moves  horizontally  to  take  the  place  of  the 
rising  air.  Thus,  a  vertically  circular  air  flow  is  established  which  stabilizes  if  an 
adiabatic  lapse  rate  is  achieved.  In  this  way,  the  large  parcels  of  air  throughout  the 
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column  are  distributed  uniformly.  Eddy  turbulenee  acts  within  the  flowing  stream  of  air, 
causing  random  swirls  and  eddies  to  arise  which  distort  the  shape  of  smaller  parcels  of  air 
within  the  stream.  Molecular  diffusion  is  the  only  process  that  does  more  than  simply 
move  pareels  of  air  around.  Taking  place  at  the  moleeular  level,  it  is  the  process  by 
which  molecules  in  the  air  (made  up  of  water  vapor  and  other  gases)  move  from  one 
parcel  of  air  to  another  due  to  differences  in  coneentration.  It  is  relatively  slow,  but 
eatalyzed  by  eddy  turbulence,  it  is  the  means  by  whieh  the  smaller  parcels  of  air  beeome 
well-mixed.  Henee,  these  three  proeesses,  spurred  by  a  variety  of  forces,  mix  the  air  to 
produce  a  vertically  homogeneous  air  column  (Livingston,  1970:  97-98). 

With  vertically  homogeneous  air  defined,  we  ean  discuss  conditions  that  cause 
departures  from  this  baseline.  Perhaps  the  most  widespread  is  heating  from  below,  whieh 
can  arise,  for  example,  when  the  ground  is  heated  in  the  morning  as  the  sun  comes  up,  or 
when  a  cool  air  mass  blows  out  over  a  warm  sea.  The  air  column,  then,  initially  in  neutral 
or  unstable  equilibrium,  eonsists  of  a  thiek  lower  region  of  cool  air  (arising,  in  the  first 
example,  from  nocturnal  cooling  of  the  surfaee,  or  in  the  seeond  example,  the  initial 
temperature  of  the  cool  air  mass),  above  whieh  the  air  is  warmer,  but  deereasing  in 
temperature  adiabatically.  The  proeess  starts  at  the  surfaee,  where  the  air  is  heated  until  it 
rises  in  patehes  through  the  eolder  air  around  it.  As  it  rises,  heavier,  eooler  air  rushes  in 
to  be  warmed  in  turn.  Thus,  convection  develops.  The  warmed  air  rises  through  the 
cooler  air  until  it  reaehes  warmer  air  above.  In  this  way,  three  regions  are  formed:  a  thin 
layer  near  the  surfaee  in  whieh  the  lapse  rate  is  greater  than  adiabatic,  but  in  which 
turbulence  is  small  because  the  air  parcels  are  just  beginning  to  accelerate;  a  thiek,  highly 
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turbulent  central  layer  through  which  the  rising  air  is  accelerating  rapidly;  and  the  upper 
region  that  is  stable.  The  process  is  illustrated  in  Figure  2.12,  where  the  numbered  curves 
indicate  successive  stages  of  the  process  either  as  the  sun  rises  (as  in  the  first  example)  or 
as  the  cool  air  mass  blows  farther  out  over  the  warm  sea  (as  in  the  second  example).  The 
intensity  of  the  effects  of  heating  from  below  is  moderated  if  there  is  cloud  cover  that 
reduces  the  surface  heating,  or  winds  that  create  so  much  surface  turbulence  they  override 
the  convective  turbulence.  Often,  these  moderating  factors  eliminate  the  homogeneous 
central  layer  altogether. 
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Figure  2.12.  Successive  temperature  distributions  due  to  heating  from  below.  Heights 
and  temperatures  are  typical,  but  can  vary  widely  (adapted  from  Kerr,  1951:  220) 

Under  typical  conditions,  however,  the  central  layer  will  be  well  mixed  with  the 
concomitant  standard  A^-profile.  The  stable  upper  layer  usually  has  a  water  vapor  lapse 
that  causes  it  to  be  superstandard,  although  occasionally  there  may  be  enough  of  a  water 
vapor  inversion  to  make  it  substandard.  The  surface  layer  is  superstandard  if  the  surface 
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is  water  or  wet  ground  and  substandard  if  it  is  dry  ground.  Over  water  especially,  a  large 
lapse  in  humidity  often  occurs  resulting  in  an  M-inversion  and  a  duct.  Heating  from 
below  is  the  foremost  cause  of  change  in  neutral  or  unstable  air. 

In  stable  air,  the  primary  change  agent  is  cooling  from  below.  This  process  occurs 
after  nightfall,  as  heat  from  the  air  diffuses  to  the  surface  directly  below  it,  or  it  can 
happen  as  a  warm  air  mass  is  blown  over  a  cool  sea.  Unlike  heating  from  below,  the  air 
is  not  vertically  turbulent,  but  becomes  more  stable.  Since  the  air  near  the  surface  is 
cooled,  it  becomes  even  more  concentrated  and  even  less  buoyant  than  it  was.  Further, 
the  vertical  movement  that  does  take  place  is  not  upwards  as  with  heating  from  below, 
but  downwards,  causing  the  temperature  and  humidity  gradients  to  be  concentrated  near 
the  surface.  Because  of  this  downward  flow  of  air,  the  overall  effects  don’t  go  as  high  as 
with  heating  from  below,  but  because  the  cool  air  is  building  up  in  the  bottom  layer, 
surface  inversions  tend  to  be  deeper  (often,  several  hundred  feet),  and,  as  time  progresses, 
they  get  continually  deeper.  Over  water  or  wet  ground,  significant  humidity  inversions 
develop  that  can  balance  out  the  temperature  inversions.  Under  these  conditions,  then, 
the  surface  layer  may  be  superstandard,  standard,  or  substandard.  Over  dry  ground,  the 
temperature  inversion  is  acting  alone,  so  the  surface  layer  is  superstandard.  Whether  over 
land  or  water,  strong  winds  and  a  small  temperature  differential  between  the  air  and  the 
surface  will  substantially  increase  the  amount  of  mixing,  reducing  the  effects  of  the 
process  on  the  refractivity  profile.  Figure  2.13  is  a  comparison  of  temperature  inversions 
under  conditions  of  strong  winds  and  small  temperature  differential,  and  light  winds  and 
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large  temperature  differential.  Thus  the  effects  of  cooling  from  below  depend  on  the 
surface  type,  the  amount  of  time  elapsed,  the  wind  speed,  and  the  temperature  differential. 

Another  factor,  which  can  affect  both  heating  from  below  and  cooling  from  below 
is  wind  shear.  Wind  shear  is  defined  as  “the  variation  with  height  of  the  horizontal 
component  of  the  wind  velocity”  (Kerr,  1951:  234).  Shear  can  be  caused  by  a  variation  in 
the  horizontal  pressure  gradient  with  height,  or,  more  usually,  by  the  influence  of  friction 
on  the  winds  at  the  surface.  Because  shear  causes  wind  speed  and  direction  to  vary  with 


(a)  (b) 

Figure  2. 13.  (a)  Cooling  from  below  with  light  winds  and  large  temperature  differential; 
(b)  Cooling  from  below  with  strong  winds  and  small  temperature  differential  (adapted 
from  Kerr,  1951:  231) 

height,  it  can  have  an  effect  on  the  refractivity  profile.  Consider,  for  instance,  warm,  dry 
land  air  that  is  blown  out  over  a  cold  sea.  Normally,  higher  winds  are  found  at  the  higher 
altitudes.  So,  in  a  given  column  of  air  over  the  sea,  the  air  higher  up  took  less  time  to  get 
into  position  and  has  been  over  the  sea  for  less  time  than  the  air  below.  For  this  reason  it 
is  warmer  than  the  air  below.  This  increases  the  thermal  stability  of  the  air  column  and 
under  typical  humidity  conditions  sets  up  a  superstandard  refracting  layer.  Another 
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example  is  that  of  a  cool  air  mass  blown  out  over  a  warm  sea.  For  the  same  reasons,  the 
column  becomes  less  stable  and  convection  is  set  up  resulting  in  a  homogeneous  layer 
that  is  close  to  standard.  Shear  also  can  alter  the  AT-profile  over  land,  especially  if  there  is 
a  marked  difference  in  surface  characteristics.  The  effects  are  seldom  as  pronounced  as 
those  near  the  coastline,  however. 

At  this  point,  at  the  risk  of  stating  the  obvious,  it  must  be  emphasized  that  the 
overwater  modifications  in  the  iV-profile  are  never  as  simple  as  the  examples  so  far  would 
suggest.  Factors  such  as  shear,  variations  in  solar  radiation,  variations  in  surface 
temperature,  and,  of  course,  lack  of  initial  homogeneity  of  the  air  column,  all  play  a  part 
in  the  final  condition  of  the  air  column.  While  the  influences  presented  above  do  play 
major  roles,  they  can  not  account  for  all  modifications  to  the  A-profile. 

In  addition  to  the  forces  acting  on  and  within  a  particular  column  of  air,  we  must 
consider  the  effects  of  the  large  air  mass  movements  within  a  geographic  region. 

One  of  these  effects  is  subsidence,  which  is  the  sinking  of  a  large  air  mass  from  a  high  to 
a  low  level.  Subsidence  is  the  result  of  dynamic  causes  and  is  generally  associated  with 
high  pressure  systems.  As  the  air  descends,  the  stable,  adiabatic,  or  unstable  nature  of  the 
air  mass  intensifies.  Because  the  subsiding  air  mass  is  usually  stable,  its  temperature  is 
usually  warmer  than  the  air  below.  Similarly,  the  air  mass  is  generally  dryer  than  the  air 
below,  since  it  descended  from  a  fairly  high  level.  These  two  factors  and  others  work  to 
produce  an  inversion  layer  between  the  sinking  air  mass  and  the  air  below  it.  This 
inversion  is  called  a  subsidence  inversion,  and  almost  always  contains  a  superstandard  M- 
gradient.  Often,  there  is  an  M  inversion,  especially  over  the  sea.  The  base  of  a 
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subsidence  inversion  is  determined  by  the  layer  of  frictional  influence  near  the  surface, 
and  can  reside  anywhere  from  one  to  four,  or  five  thousand  feet  up.  The  thickness  of  a 
subsidence  inversion  can  be  several  hundred  feet  or  a  several  thousand  feet.  Obviously, 
an  inversion  this  large,  under  the  right  conditions,  will  severely  impact  radio  propagation. 
Ducting  is  extremely  common  near  the  great  high  pressure  systems  that  occur  over  the 
oceans  around  the  world  at  about  30  degrees  north  and  south  latitude.  Near  the  southwest 
coast  of  California  ducts  occur  on  an  average  of  40  percent  of  the  time,  and  along  the 
coast  of  Japan,  about  10  percent  of  the  time  (Patterson,  1994:  16).  Additionally, 
subsidence  inversions  can  also  occur  on  the  lee  side  of  mountains. 

When  large  air  masses  collide  horizontally,  they  form  a  front.  An  air  mass,  by 
definition,  is  approximately  uniform  with  regard  to  temperature  and  humidity.  When  two 
air  masses  traveling  horizontally,  come  into  contact  with  each  other,  the  warmer  air  mass 
slides  up  over  the  cooler  one.  The  boundary,  or  front,  between  the  masses  is  always  close 
to  horizontal.  Though  its  slope  increases  with  wind  speed,  increasing  latitude,  and 
temperature  similarity,  it  generally  remains  somewhere  between  0  and  1/25.  Because  the 
upper  air  mass  is  warmer,  there  is  then  always  a  stable  layer  at  the  front,  and  often  there  is 
a  temperature  inversion.  Ordinarily,  the  water  vapor  concentration  increases  with  height 
in  this  layer,  but  not  always.  For  this  reason,  the  front  may  define  a  superstandard  layer,  a 
standard  layer,  or  a  substandard  layer.  Frontal  inversions  do  not  lend  themselves  to 
generalities.  The  passing  of  a  front  must  also  be  regarded  with  great  care  because  it 
nearly  always  precedes  a  potentially  radical  change  in  weather  conditions  across  the 
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landscape  and  over  time,  which  in  turn  causes  the  M-profile  to  vary  with  distance  and 
time. 

A  frontal  boundary  is  but  one  of  many  elements  that  will  affect  the  horizontal 
variation  of  the  refractivity  profile.  Though  this  thesis  is  primarily  concerned  with  long 
term  averages  and  climatological  predictions  of  the  vertical  variations,  the  radar  engineer 
must  take  into  account  all  the  local  conditions,  either  current  or  predicted,  to  wisely 
understand  how  the  climatological  averages  may  vary.  Consideration  of  the  horizontal 
gradient  is  inescapable.  Though  coastlines  probably  produce  the  most  outstanding 
horizontal  aberrations,  similar,  if  less  intense,  phenomena  may  be  found  wherever  surface 
characteristics  or  weather  systems  change. 

None  of  the  forces  described  above  acts  independently.  There  are  always  larger 
influences  driving  theses  forces,  and  while  we  can  rarely  predict  exactly  what  weather 
conditions  will  occur  and  exactly  what  effect  they  will  have  on  the  refractivity  profile,  we 
can,  by  studying  climatic  trends,  determine  how  geography,  topology,  diurnal  cycles,  and 
seasonal  cycles  have  traditionally  impacted  the  local  weather  patterns  at  some  time  of 
year  and  time  of  day.  Here  are  some  of  the  major  climatic  effects  followed,  in  the  next 
section,  by  a  survey  of  worldwide  climate  types. 

Inarguably  the  most  regular  and  predictable  climatic  effect  is  the  diurnal  cycle, 
that  is,  the  regular  procession  of  day  to  night  and  back  to  day.  Assuming  no  irregular 
influences,  the  cycle  will  proceed  something  like  this.  At  midday,  the  atmosphere  has 
been  fully  heated  with  all  the  accompanying  convective  mixing,  so  it  will  be 
homogeneous.  As  sunset  approaches,  the  ground  begins  to  cool  as  the  sun’s  rays  become 
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less  and  less  direct.  Cooling  from  below  starts  and  stratification  begins  to  develop  near 
the  surface,  but  the  Af-profile  is  probably  affected  only  very  little  until  after  sunset. 
During  the  night,  the  surface  cools,  radiating  heat.  Much  of  this  heat  escapes  through  the 
atmosphere,  but  some  of  it  is  absorbed  and  reradiated  by  atmospheric  gases,  especially 
water  vapor.  For  this  reason,  the  presence  of  clouds  reduces  the  amount  of  surface 
radiation  considerably.  At  this  point,  cooling  from  below  over  water  continues  pretty 
much  as  described  in  that  section.  Cooling  of  the  air  over  land,  however,  is  more 
complicated  both  with  regard  to  the  rate  and  amount  of  cooling,  and  to  how  the  air  above 
the  surface  is  affected.  First,  the  temperature  of  land  varies  greatly,  both  with  time  as  the 
night  wears  on,  and  topology,  if  the  surface  characteristics  vary.  Temperature  variations 
in  water  are  much  less  intense.  In  particular,  the  temperature  and  amount  of  radiation 
from  land  depend  upon  conductivity  of  the  soil,  which  determines  how  fast  the  radiating 
heat  is  replaced  from  below.  For  example,  dry  sand  is  about  three  times  as  conductive  as 
dry  soil,  and  wet  soil  is  almost  three  times  as  conductive  as  dry  sand.  Furthermore,  land, 
unlike  water,  does  not  uniquely  determine  the  amount  of  humidity  at  the  surface,  unless 
the  ground  is  wet.  Immediately  over  water,  the  air  is  saturated;  but  over  land  humidity 
can  vary  greatly.  On  the  other  hand,  if  the  land  is  completely  dry,  there  will  be  no 
humidity  gradient  at  all,  and  the  M-profile  will  depend  solely  upon  temperature.  Because 
of  the  variability  of  these  factors,  a  variety  of  M-profiles  may  occur  at  night.  Both 
superstandard  and  substandard  anomalies  may  appear  both  at  the  surface  and  aloft.  The 
superstandard  behavior  may  produce  ducts  that  grow  in  the  evening  just  after  sunset, 
intensify  during  the  night,  and  dissipate  with  the  sunrise.  The  most  important 
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consideration,  however,  is  that  diurnal  anomalies  primarily  occur  at  night.  In  the  morning, 
the  sun  again  warms  the  earth,  heating  from  below  begins,  and  through  eonveetion  and 
turbulent  mixing  the  atmosphere  again  becomes  homogeneous. 

The  wild  anomalous  changes  in  the  M-profile  that  occur  during  the  diurnal  cycle 
are  generally  more  intense  during  the  summer  than  the  winter.  There  are  three  primary 
reasons.  First,  the  ground  is  typically  drier  in  the  summer  than  in  the  winter.  Sinee  dry 
ground  has  a  lower  conductivity,  the  surface  temperature  changes  in  the  summer  are  more 
dramatie  than  in  the  winter,  producing  stronger  temperature  inversions.  The  second 
reason  is  that  the  water  vapor  content  is  usually  higher  in  the  warm  summer  air,  resulting 
in  larger  humidity  gradients.  Finally,  the  cloudiness  and  strong  winds  that  moderate  the 
diurnal  cycle,  are  more  likely  to  occur  in  winter  than  in  summer,  at  least  in  the  temperate 
latitudes. 

Another  climatologieal  phenomenon  that  is  widespread  enough  to  be  eovered 
separately  is  the  sea  breeze.  Since  many  radars  operate  along  the  coastline,  what 
happens  there  is  particularly  important  to  understand.  Sea  breezes  arise  during  the  day 
when  warm  land  air  rises  and  is  subsequently  replaced  by  cool  sea  air  rushing  in  to  fill  the 
void.  In  this  way  a  cireulation  from  sea  to  land  is  set  up.  When  the  land  air  has  traveled 
some  distance  out  to  sea  it  subsides,  often  setting  up  a  subsidenee  inversion  depending 
upon  the  humidity  distribution.  This  inversion  may  result  in  dueting.  Similarly,  at  night 
the  land  eools  faster  than  the  sea  eausing  the  warm  sea  air  to  rise  over  the  cooler  land  air 
creating  a  cireulation  in  the  opposite  direction  called  a  land  breeze.  Here  again,  dueting 
can  occur  because  of  subsidenee. 
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Finally,  it  must  be  noted  that  extensive  snow  cover,  fog,  and  cloudiness  will  all 
act  to  moderate  to  some  extent  the  effects  described  in  this  section.  Snow  cover  reduces 
both  solar  heating  of  the  ground  during  the  day  and  radiation  from  the  ground  at  night, 
especially  in  the  high  latitudes.  Cloud  cover  will  have  a  similar  effect  because  it 
moderates  both  solar  heating  and  nocturnal  radiation.  When  fog  forms,  water  vapor  in 
the  air  changes  to  liquid,  reducing  the  vapor  pressure.  The  humidity  lapse  that  results 
tends  to  counteract  any  temperature  inversions  that  exist,  preventing  extreme 
superstandard  behavior  (Bean  and  Dutton,  1966:  135). 

Climate  Survey 

While  the  meteorological  effects  described  in  the  last  section  determine  the 
anomalous  characteristics  of  the  M-profile,  the  value  of  N  at  the  surface  and  the  initial  N- 
gradient  apart  from  anomalous  influences  are  determined  primarily  by  climate,  that  is 
geography  and  time  of  year.  Based  on  over  two  million  weather  observations  taken  at  45 
weather  stations  across  the  continental  United  States,  the  average  value  of  refractivity  at 
the  surface,  Ns,  is  313,  according  to  Bean  and  Thayer  (1959).  Using  Bean  and  Thayer’s 
CRPL  standard  atmosphere  (2-10),  we  can  calculate  the  initial  gradient,  AN,  to  be 
approximately  -41.94.  These  numbers  represent  the  average.  However,  seasonal  changes 
and  the  vast  array  of  topographies  across  the  United  States  give  us  surface  values  that 
vary  widely  from  these.  When  the  scope  of  observation  is  expanded  worldwide,  the 
variations  become  even  more  expansive. 
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Bean  and  Dutton  (1966:  89-109)  conducted  a  study  of  how  surface  refract! vity 
varies  with  climate  using  long  terms  means  of  temperature,  pressure,  and  humidity  given 
in  the  United  Nations  monthly  publication  Climate  Data  for  the  World.  They  examined 
five  years  of  records  for  each  of  306  stations  worldwide  from  the  period  1949-1958.  This 
section  contains  a  summary  of  their  most  significant  findings. 

Most  of  Bean  and  Dutton’s  conclusions  concerned  the  annual  range  of  variation  of 
surface  values  of  N.  Ns  varies  with  the  diurnal  cycle,  with  season,  and,  of  course,  with 
topology.  Table  2.  2  contains  the  annual  means  and  ranges  of  Ns  for  six  different  climate 
types  which  are  defined  and  typified.  From  the  data  presented  there,  the 


Table  2.  2.  Characteristics  of  climatic  types  (Bean  and  Dutton,  1966:  103) 


Type 

Location 

Annual 
mean  Ns 
(N  units) 

Annual 
range  of  Ns 
(N  units) 

Characteristics 

I.  Midlatitude- 
coastal 

Near  the  sea  or  in  lowlands  on  lakes  and 
rivers,  in  latitude  belts  between  20°  and 
50°. 

300  to  350 

30  to  60 

Generally  subtropical  with 
marine  or  modified  marine 
climate. 

II.  Subtropical- 
Savanna 

Lowland  stations  between  30°N  and 

25 °S,  rarely  far  from  the  ocean. 

350  to  400 

Definite  rainy  and  dry 
seasons,  typical  of 

Savanna  climate. 

III.  Monsoon- 
Sudan 

Monsoon— generally  between  20°  and 
40°N,  Sudan-across  central  Africa 
from  10°  to  20°N 

280  to  400 

60  to  100 

Seasonal  extremes  of 
rainfall  and  temperature. 

IV.  Semiarid- 
Mountain 

In  desert  and  high  steppe  regions  as 
well  as  mountainous  regions  above 

3,000  ft. 

240  to  300 

0  to  60 

Year-round  dry  climate. 

V.  Continental- 
Polar 

In  middle  latitudes  and  polar  regions. 
(Mediterranean  climates  are  included 
because  of  the  low  range  resulting  from 
characteristic  dry  summers). 

300  to  340 

0  to  30 

Moderate  or  low  annual 
mean  temperatures. 

VI.  Isothermal- 
equatorial 

Tropical  stations  at  low  elevations 
between  20°N  and  20°S,  almost 
exclusively  along  seacoasts  or  on 
islands 

340  to  400 

Oto  30 

Monotonous  rainy 
climates. 

following  observations  can  be  made:  Generally,  Ns  decreases  with  increasing  latitude. 
Surface  refractivity  values  in  purely  maritime  regions  such  as  the  west  coasts  of  North 
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America  and  Europe  are  low  and  undergo  relatively  little  variation  due  to  the  constant 
onshore  advection  of  cool,  moist  sea  air.  Continental  areas  like  central  North  America 
generally  have  surface  refractivity  lower  than  the  adjacent  maritime  areas,  but  with  more 
variation.  The  variation  comes  from  the  changes  in  radiative  heating  from  day  to  night, 
and  summer  to  winter,  which  are  not  modified  by  air  from  the  ocean.  Regions  similar  to 
the  east  cost  of  the  United  States  have  a  combined  maritime  and  continental  influence,  in 
which  a  variety  of  air  mass  types  interact  to  produce  a  moderate  amount  of  variation 
throughout  the  year.  Mountainous  places  generally  have  lower  values  of  Ns  because  of 
the  fundamental  changes  in  temperature,  pressure,  and  humidity  with  altitude.  The 
refractivity  variation  is  tj^ically  moderate  in  the  mountains.  Larger  ranges  are  found  in 
places  like  Australia,  the  African  plateau  near  the  Cameroons,  and  in  the  Great  Basin  of 
the  southwestern  United  States  where  there  are  wide  ranging  temperature  variations 
throughout  the  year,  but  relatively  little  humidity  to  act  as  a  moderating  factor.  Some  of 
the  largest  ranges  in  the  world  are  found  in  the  African  Sudan  and  in  areas  influenced  by 
the  unique  characteristics  of  the  Indian  monsoon. 

Of  course,  this  has  been  an  extremely  broad  look  at  world  climate  types.  Studies 
of  local  patterns  in  a  given  region  would,  of  course,  provide  more  insight  into  the 
refractive  effects  playing  in  that  area.  Indeed,  Bean  and  Dutton  (1966:  1 10-131)  include  a 
fairly  detailed  overview  of  refractivity  in  the  continental  United  States  that  may  be  of 
interest  to  the  radar  engineer  involved  in  refractivity  prediction. 
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With  a  detailed  knowledge  of  worldwide  climate,  a  solid  grasp  of  the  various 
forces  that  can  alter  the  M-profile,  and  a  thorough  understanding  of  the  mathematics  of 
N-refractivity  and  raytracing,  one  may  construct  a  model  that  can  predict,  with  some 
accuracy,  how  much  the  atmosphere  will  bend  an  electromagnetic  radar  wave  launched 
into  the  atmosphere  at  a  given  elevation  angle.  Accuracy,  here,  is  relative  and  must  be 
compared  to  what  has  gone  before.  Previous  radar  range  performance  models  have  relied 
on  four-thirds  earth,  or  at  best,  a  standard  refractive  profile  based  on  the  United  States 
average  surface  refractivity  value,  Ns=313.  Climatology  and  raytracing  can  put  us  closer. 
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in.  Methodology 


Construction  of  a  refractivity  profile  depends  on  how  much  knowledge  there  is  of 
the  prevailing  atmospheric  conditions.  Over  the  past  150  years,  sailors,  scientists  and 
engineers  have  collected  a  wealth  of  climatological  data  fit  for  making  educated  guesses 
about  what  atmospheric  conditions  will  be  like  at  a  given  place  and  time.  A  profile  built 
up  from  this  type  of  data  is  all  that  is  needed  to  compute  a  close  approximation  to  the 
most  direct  path  from  a  given  radar  to  a  given  target  on  a  given  day  in  a  particular 
location. 

Specifically,  geometric  optics  (raytracing)  is  used  to  calculate  the  path  of  a  single 
ray  of  energy  leaving  the  radar  at  some  takeoff  (initial  elevation)  angle.  The  only 
problem  is  determining  which  takeoff  angle  will  put  the  endpoint  of  the  ray  on  the  target. 
Since  the  raytrace  is  not  a  closed-form  solution,  there  is  no  way  of  working  backward  to 
the  answer.  Hence,  a  trial-and-error  iterative  approach  is  employed:  First,  an  initial 
takeoff  angle  is  estimated  and  the  ray  is  traced  through  the  atmosphere  until  it  reaches  the 
target  height.  Next,  the  subtense  of  the  trace  endpoint  (i.e.  the  angle  defined  by  the  radar, 
the  earth’s  center  and  the  endpoint  of  the  trace)  is  compared  with  the  subtense  of  the 
target,  and  a  new  takeoff  angle  is  estimated  based  on  the  magnitude  and  direction  of  the 
difference.  The  process  is  repeated  until  the  trace  ends  at  the  target  subtense  (Figure  3.1). 
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Figure  3.1  Initial  takeoff  angle,  ai,  estimated;  trace  1  goes  beyond  target.  New  angle, 
az,  estimated;  trace  2  undershoots.  Endpoint  of  trace  3,  using  angle,  as,  finds  target. 

Finally,  the  initial  elevation  and  the  path  length  (or  apparent  range)  are  used  to  calculate 
the  apparent  height,  which,  when  compared  with  the  actual  height,  yields  the  height  error. 
Similarly,  the  trace  length  (apparent  range)  compared  with  the  actual  geometric  (or 
straight  line)  range  provides  the  range  error. 

This  chapter  details  the  above  process.  First,  the  database  is  described.  Following 
that  is  a  step-by-step  discussion  of  how  the  model  developed  for  this  project  works.  Note 
that  throughout  the  discussion,  the  symbols  used  for  the  variables  involved  match  as 
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closely  as  possible  those  used  in  the  MATLAB®  code  itself.  Also  note  that  the  important 
capabilities  and  limitations  of  the  model  are  specifically  discussed  in  Chapter  5,  Findings 
and  Conclusions.  Also,  a  complete  listing  of  the  MATLAB®  source  code  is  included  in 
Appendix  A. 

Climatology  Database 

The  fuel  for  this  climatology-based  radar  propagation  prediction  model  is  a  the  Historical 
Electromagnetic  Propagation  Condition  (HEPC)  Database.  Compiled  by  the  U.S.  Navy, 
the  HEPC  database  contains  statistical  climatological  data  for  921  radiosonde  stations 
distributed  worldwide  and  is  an  integral  part  of  the  Navy’s  Advanced  Refractive  Effects 
Prediction  System  (AREPS)  propagation  model.  The  statistics  it  contains  are  derived 
from  two  meteorological  databases:  GTE  Sylvania’s  Radiosonde  Data  Analysis  II,  and  the 
National  Climatic  Data  Center’s  Duct63  database.  The  former  is  a  “large  scale  analysis 
of  approximately  three  million  worldwide  radiosonde  soundings  from  1966  to  1969  and 
1973  to  1974.”  The  latter  consists  of  mostly  marine  surface  observations  spanning  15 
years.  These  observations  were  taken  from  ship  logs,  ship  weather  reporting  forms, 
published  ship  observations,  automatic  buoys,  teletype  reports,  and  card  decks  purchased 
from  foreign  meteorological  services  (Patterson,  1987: 1-3,9). 

The  HEPC  database  is  public  domain,  available  on  the  Space  and  Naval  Warfare 
Systems  Center  webpage,  and  is  encoded  in  .dbf  format  which,  with  the  exception  of  a 
header  and  some  delimiters,  is  ASCII  text. 
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Each  record  in  the  database  contains  the  statistics  for  a  single  radiosonde  station. 
The  first  seven  fields  in  a  record  are  station  information.  The  next  twelve  each  contain 
monthly  mean  refractivity  parameters.  The  fields  are  defined  as  follows  (Figure  3.2): 

MSQ  -  Marsden  Square  containing  the  station.  A  Marsden  Square  is  a  10°  lat.  by 
10°  long,  portion  of  the  Earth’s  surface 
WMO  -  World  Meteorological  Organization  station  number 
Name  -  station  name 

Inland/Coast  -  Indicator  (L=inland,  C=coast) 

Latitude  -  sign  is  negative  for  south,  positive  for  north 
Longitude  -  sign  is  negative  for  east,  positive  for  west 
Station  Elevation  -  above  mean  sea  level  (meters) 

Monthly  Data  (one  field  for  each  month)  -  Each  contains  the  following: 

Number  of  1200Z  observations 

Number  of  OOOOZ  observations 

Surface  N-unit  value 

Surface  to  1000  meter  M-unit  gradient 

Surface-based  duct  M-unit  gradient 

Sfc-based  duct  optimum  height  (layer  base)  (meters) 

Surface-based  duct  thickness  (meters) 

Surface-based  duct  M-unit  deficit 
Surface-based  duct  trapping  frequency  (MHz) 

Sfc-based  duct  1200Z  percent  occur 
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Sfc-based  duct  OOOOZ  percent  occur 


Elevated  duct  M-unit  gradient 

Elevated  duct  optimum  height  (layer  base)  (meters) 

Elevated  duct  thickness  (meters) 

Elevated  duct  M-unit  deficit 
Elevated  duct  trapping  frequency  (MHz) 

Elevated  duct  1200Z  percent  occur 
Elevated  duct  OOOOZ  percent  occur 
Probability  of  >  1  elevated  duct  (%xl00) 

Probability  of  sfc-based  and  elevated  ducts  (%xl00)  (Patterson,  1993) 
Statistical  data  on  surface  conditions  and  the  significant  ducting  parameters  is 
listed  this  way  for  every  month  of  the  year  for  virtually  every  place  in  the  world  inhabited 
by  human  beings  (and  some  that  are  not).  Sadly,  the  database  contains  only  mean  values 
for  each  parameter.  There  is  no  variation  data.  This  omission  certainly  obscures  the  true 
climatological  picture,  but  as  no  other  data  set  is  as  complete  and  usable,  the  partial  view 
must  be  accepted. 

Each  data-element  described  above  is  used  in  CLIMAREF  at  some  point.  Indeed, 
the  database  is  consulted  throughout  the  model,  from  the  initial  stages  of  determining  the 
most  appropriate  radiosonde  station  to  use,  to  the  construction  of  the  refractivity  profile 
before  the  final  raytracing. 
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-  Each  record  consists  of  a  header  and  12  month-fields: 

(Header) 


MSQ 

WMO 

Name 

Inland/Coast  indicator 
(L==inland,  C=coastal) 

Latitude 

Longitude 

Elevation 

I  062 

48354 

Udom  Thani,  Thailand 

L 

17.37 

-102.80 

178 

(Jan  month  field  -  1  exists  for  each  month  of  the  year) 
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Figure  3.2  Fields  contained  in  one  record  of  the  HEPC  database 


CLIMAREF  Detailed  Description 

CLIMAREF,  the  model  constructed  for  this  thesis,  was  written  in  MATLAB®, 
Version  5.2.0.3084.  There  are  nine  modules,  or  in  MATLAB®  terminology,  “functions.” 
(see  Appendix  A).  The  first  module  calls  the  other  eight  and  performs  various  display 
and  plotting  functions.  This  module  (and  possibly  module  2,  the  input  routine)  may  be 
tailored  or  replaced  in  order  to  integrate  CLIMAREF  into  a  larger  application  without 
having  to  modify  the  other  eight  modules. 

The  remainder  of  the  chapter  contains  a  theoretical  description  of  how 
CLIMAREF  works.  While  the  operation  of  every  module  is  covered,  there  is  no  attempt 
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made  beyond  this  point  to  delineate  where  one  module  leaves  off  and  another  begins. 

The  MATLAB®  routines  are  internally  documented  fairly  well,  and  are  included  in  the 
Appendix  A. 

Initial  Data  Entry.  The  input  data  the  model  requires  are  month  of  operation, 
MON,  whether  operation  will  be  during  the  day  or  night,  DYNT,  the  latitude  and 
longitude  of  the  radar,  RLAT,  RLONG,  the  elevation  of  the  radar  above  ground  level, 

REL,  the  direction  the  radar  is  pointing  in  degrees,  RAZ,  the  height  of  the  target  above 
ground  level,  THT,  and  the  range  of  the  target,  TRG. 

Two  important  qualifications  regarding  this  data  must  be  made.  First,  all  heights 
and  ranges  must  be  converted  to  kilometers  before  they  are  used  in  the  program.  The  data 
may  be  entered  in  feet  or  nautical  miles,  but  it  must  be  converted.  Second,  all  heights  are 
measured  firom  ground  level,  defined  as  the  elevation  above  sea  level  of  the  radiosonde 
station  whose  data  is  used  to  construct  the  N-profile.  It  is  assumed  the  radar  will  be  near 
enough  to  the  station  that  surface  elevation  of  the  two  locations  will  be  fairly  similar.  If 
not  the,  station  selected  may  not  be  the  best  station  for  the  purpose  even  if  it  is  the 
nearest.  More  will  be  discussed  on  this  subject  later. 

Choosing  station(s)  to  use  for  calculation.  Given  the  radar  latitude  and  longitude, 
CLIMAREF  finds  the  ten  nearest  radiosonde  stations  to  the  radar  site.  The  user,  given 
the  list,  is  then  prompted  to  choose  the  best  station  for  the  purpose.  To  identify  the 
nearest  stations,  CLIMAREF  takes  advantage  of  the  Marsden  Square  method  of 
organizing  the  stations. 
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Each  radiosonde  station  in  the  HEPC  database  is  labeled  with  a  WMO  number 
and  a  Marsden  Square  number.  The  WMO  number  is  simply  for  identification.  The 
MSQ  number,  however,  locates  that  station  in  a  square,  so  to  speak,  bordered  by  parallels 
and  meridians  evenly  divisible  by  ten,  e.g.  Dayton  is  in  MSQ  1 17,  bordered  by  the  30th 
and  40th  parallels  and  the  80th  and  90th  west  meridians.  Further,  there  are  limits  to  MSQ 
coverage.  There  are  no  MSQs  defined  north  of  80°  N  or  south  of  70°  S.  Within  that  area, 
the  radiosonde  stations  are  most  concentrated  near  populated  areas,  espeeially  in  the  more 
developed  countries. 

Using  the  MSQ  numbering  rules  (Patterson,  1987:  10)  and  a  look-up  table  of  the 
MSQs  adjacent  to,  and  west  of,  the  prime  meridian,  CLIMAREF  determines  the  number 
of  the  MSQ  where  the  radar  is  located  and  the  number  of  each  of  the  surrounding  MSQs 
(Figure  3.3).  Next,  using  pointers  collected  from  a  reference  file  called  MSQLIST, 
CLIMAREF  identifies  all  the  stations  located  within  those  nine  MSQs  and  gathers 
location  data  for  each  station.  MSQLIST  was  created  for  this  projeet  as  a  quiek-reference 
companion  file  to  the  main  database. 

To  determine  which  stations  of  those  within  the  nine  MSQs  are  the  nearest  to  the 
radar,  a  spherical  coordinate  system  is  assumed  at  the  origin  of  the  earth  and  the  position 
of  each  of  the  stations  is  determined  based  on  the  latitude  complement,  the  longitude,  and 
the  local  earth  radius  of  curvature  (Figure  3.4).  The  latter  is  ealculated  (Abel,  1982: 
18-19): 

_ _ 

^  5(l-hCcOS^(|))'^^(l-|-CcOS^(t)COS^0)  ’ 
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where. 


A  =  semimajor  (equatorial)  radius  (6378.139  km) 

B  =  semiminor  (polar)  radius  (6356.750  km) 

C  =  /  5^  - 1 

(])  =  latitude 
0  =  radar  azimuth 

The  local  radius  of  curvature  is  used  because  the  earth  more  closely  approximates  an 
oblate  spheroid  than  a  true  sphere.  Of  course,  the  curve  of  the  local  earth  surface  is  the 
important  factor  in  the  path  calculations,  not  the  actual  radius. 


Figure  3.3  Radar  (star)  and  surrounding  radiosonde  stations  (circles)  within  the 
nine-MSQ  group 
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Then,  the  spherical  coordinates  are  converted  to  three-dimensional  rectangular  using  the 
standard  formulas: 


x  =  RCsin{RLONG)cos{90-RLAT)  ,  (3-2) 

y  =  RCsmiRLONG)cos{90-RLAT),  (3-3) 


z^RCcosiRLONG)  .  (3-4) 

If  we  follow  the  same  procedure  to  determine  the  coordinates  of  the  radar,  we  can 
compare  each  station  vector  to  the  radar  vector  to  get  the  subtense  angle. 


cos  P/J5  -  ^^2 


(3-5) 


where, 

Xs,ys,Zs  =  station  coordinates 
XR,y^,z^  =  radar  coordinates 

=  subtense  between  radar  and  station 


Figure  3.4  Calculating  the  distance  between  the  radar  and  a  nearby  radiosonde  station 
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and  then  the  actual  surface  distance  between  them  is: 


distj^s  =  X  RC.  (3-6) 

There  is  some  small  inaccuracy  in  this  calculation  since  the  earth  only  approximates  an 
oblate  spheroid,  but  for  these  purposes  it  is  quite  acceptable. 

Once  the  distance  to  every  station  within  the  nine-MSQ  region  is  calculated  the 
list  of  stations  is  sorted  by  distance  and  truncated  after  ten  stations.  Ten  is  an  arbitrary 
cutoff  with  no  special  significance  so  five  or  some  other  number  might  work  just  as  well. 
The  important  thing  is  that  there  be  enough  stations  listed  to  allow  the  user  some 
flexibility  in  his  choice. 

The  final  step  in  the  station  selection  process  is  to  display  these  ten  nearest 
stations  along  with  the  elevation,  and  distance  from  the  radar  for  each,  allowing  the  user 
to  choose  one  or  more  stations  with  which  to  approximate  the  atmospheric  conditions 
local  to  his  radar.  Normally,  the  user  will  choose  the  nearest  station.  If,  however,  the 
radar  is  situated  in  the  midst  of  two  or  three  stations  which  all  contain  relevant  data,  the 
user  may  seleet  both  or  all  of  them  and  let  the  model  interpolate  the  data  (though  at 
present  the  model  will  not  interpolate  ducting  data  —  only  surface  parameters).  In  many 
cases  this  will  be  the  most  accurate  option.  Care  must  be  taken,  however,  whether 
choosing  a  single  station  or  interpolating,  that  the  climate  of  the  radar  and  the  climate  of 
the  station(s)  are  similar.  For  instance,  consider  a  radar  at  1000  feet  above  mean  sea  level 
(MSL);  the  nearest  station  is  150  km  away  at  500  feet  MSL;  and  the  next  station  is  180 
km  away  at  900  feet  MSL  on  the  same  plateau  as  the  radar.  In  this  ease  it  is  a  bad 
decision  to  choose  the  nearer  station  because  of  the  altitude  difference.  The  second 
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station  will  likely  have  a  climate  more  closely  approximating  that  of  the  radar.  This  same 
reasoning  will  apply  when  choosing  two  or  three  stations  for  interpolation.  Interpolating, 
for  example,  between  an  off-coast  (fixed  ship)  station,  and  an  inland  site  is  hardly  a  wise 
idea.  Because  of  the  abstract  nature  of  this  final  decision,  it  is  left  to  the  user  to  make. 

Constructing  the  Refractivity  Profile.  When  the  decision  is  made,  and  the  user 
selects  the  indices  of  one,  two,  or  three  stations,  CLIMAREF  will  retrieve  the  data  from 
the  database,  interpolate  if  necessary,  and  calculate  the  refractivity  profile,  N(h),  from  the 
surface  to  the  highest  significant  altitude. 

For  each  station  selected,  CLIMAREF,  retrieves  the  surface  refractivity,  NS,  and 
the  modified  refractivity  gradient  up  to  1000  meters,  AMIK.  This  quantity  is  converted  to 
N-units  (Patterson,  1987:  14), 

10^ 

iSN\K  =  m\K- - =  AMlii:-156,  (3-7) 

where  Og  is  the  mean  radius  of  the  earth,  and  used  to  compute  the  refractivity  at  1000 
meters: 

N\K=NS  +  NN\K.  (3-8) 

These  two  parameters,  NS,  and  NIK,  are  the  basis  for  construction  of  the  refractivity 
profile  and  must  be  obtained  for  each  of  the  one,  two,  or  three  stations  selected. 

Additionally,  a  surface  duct  or  an  elevated  duct  may  be  included  in  the  refractivity 
profile.  If  the  decision  is  made  not  to  interpolate  and  a  single  station  is  selected,  the 
ducting  statistics  are  retrieved  from  the  database  and  displayed  along  with  a  query 
requiring  the  user  to  select  “Surface  duct,  Elevated  duct,  or  No  duct.”  The  user  may 
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make  his  selection  based  on  the  displayed  “Percent  of  time  duct  occurs”  and  the  heights 
and  thicknesses  associated  with  the  duct.  Remember  again,  that  these  parameters  are  only 
mean  values.  Anomalous  propagation  effects  are  highly  unpredictable  and  the  mean 
parameters  are  certainly  nothing  to  be  counted  on.  Nonetheless,  some  users  may  want  to 
include  ducting  to  get  an  idea  of  the  effects  a  duct  would  have  if  present.  The  extent  to 
which  ducting  is  successfully  modeled  is  discussed  in  the  next  two  chapters. 

If  the  user  chooses  more  than  one  station  so  as  to  arrive  at  an  interpolated  profile, 
he  will  forego  ducting  by  default.  There  are  two  interpolation  routines,  one  for  the  two- 
station  scenario,  and  one  for  the  three-station  scenario.  For  the  purposes  of  this  paper,  the 
former  will  be  called  linear-interpolation,  and  the  latter,  planar-interpolation.  Both 
compute  new  values  of  NS,  and  NIK.  Both  of  these  algorithms  were  custom-built,  so  to 
speak,  for  this  project,  and  will  be  derived  in  some  length  here. 

The  linear-interpolation  algorithm  is  simple  as  interpolations  go  (Figure  3.5a). 
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Figure  3.5  (a)  Linear  and  (b)  Planar  interpolation  methods 
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First,  the  distance  between  the  two  stations  is  calculated  in  the  same  way  the  distances 
between  the  radar  and  the  stations  were  calculated.  Then,  a  triangle  is  constructed  using 
the  distances  between  the  three  sites  (i.e.  radar,  station  1,  station  2):  The  radar,  point  r,  is 
placed  on  the  Cartesian  x-y  plane  at  (0,0),  and  station  1,  point  si,  at  (<i;,0),  dj  being  the 
distance  between  the  radar  and  station  1.  The  coordinates  of  station  2,  point  s2,  are 
calculated  using  the  other  two  distances  (Figure  3.6): 

xl+yl=dl  (3-9) 


and 


hence. 


(<^1  •^)  +  >'2  —  ^12 


dl-di2  +dl 


(3-10) 


2^1 


3^2 


=  ^|d■ 


2  '^2  • 


(3-11) 

(3-12) 


Note  it  doesn’t  matter  whether  station  2  is  placed  in  the  positive  or  the  negative 
y-halfplane. 


Figure  3.6  Geometry  for  x-y  placement  during  linear  interpolation 
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The  z-coordinates  of  each  station  will  be  NS,  eind  then  NIK,  depending  upon 
which  value  is  being  interpolated.  First,  however,  these  values  must  be  normalized  to  sea 
level  and  sea  level  +  1km,  respectively,  since  the  elevations  of  the  two  radiosonde 
stations  will  probably  be  different.  In  other  words,  we  must  find  a  common  height  at 
which  to  interpolate  the  refractivity  values.  This  is  done  using  the  bi-exponential 
refractivity  model  (see  eq.  2-8, 2-9): 


N, 


(2-9) 


\km 


N  =  Nsexp{-c^(h-hs)},  (2-8) 

Nmsl  =  exp{-c,[fl,  -  (a,  +  SEL)] }  =  Ns  exp{-c^[-5EL] } ,  (3-13) 

where  h  is  the  height  at  which  we  want  to  find  N,  hs  is  the  surface  height,  is  the  mean 
earth’s  radius,  and  SEL  is  the  station  elevation.  Thus  we  have  two  points,  Sin  and  S2n,  in 
three-dimensional  space  representing  the  locations  and  IV-values  of  the  two  stations,  and  a 
set  of  x-y  coordinates,  r,  representing  the  location  of  the  radar  site. 

To  perform  the  interpolation  itself,  a  line  is  constructed  between  sm  and  S2n: 


x-Xi  ^  y-y^  ^  z-z^ 
a  b  c 


(3-14) 


where. 


a  =  X2-x^  b  =  y2-yi  c  =  Z2-Zi  ■  (3-15) 

Next,  we  consider  the  projection  of  the  line  onto  the  x-y  plane,  and  find  the  slope  of  the 
projection  to  be  m  =  bla.  We  want  to  find  the  value  of  N  (surface  or  IK)  at  the  point,  p, 
on  the  line  nearest  to  the  radar  (Figure  3.7).  So,  we  construct  a  perpendicular  line 
(m  =  -alb)  from  line  S]-S2  to  the  radar  location,  point  r. 
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Figure  3.7  Construction  of  perpendicular  during  linear  interpolation 
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-x,  +  -x^-y^  +  y^ 
b  a _ 

b  a 
a  b 


(3-17) 


As  it  turns  out,  we  need  not  calculate  yp.  To  get  our  interpolated  value  we  simply 
calculate  Zp  ,  the  height  value  of  line  S1-S2  at  point  p: 


Z.-Zl  X  -x^ 


a 


interp  * 


(3-18) 


Finally,  the  radar  surface  elevation  is  interpolated  from  the  station  heights  and  the 
interpolated  values  of  NS  and  NIK  (at  sea  level)  are  de-normalized  to  the  interpolated 
elevation  using  equations  3-13  and  3-14. 

The  three-station,  or  planar,  interpolation  is  similar  to  the  two-station  case  except 
that  a  plane  is  constructed  between  the  stations  instead  of  a  line  (Figure  3.5b).  First  a 
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triangle  is  constructed  as  in  the  two-station  case  between  points  r,  Si,  and  S2  (Figure  3.8). 
Then  xs  and  ys  are  computed  the  same  as  X2  and  y2  (equations  3-1 1,  3-12). 


83(^3.  ya) 


Figure  3.8  Geometry  for  x-y  placement  during  planar  interpolation 


Again,  y2  is  considered  positive  as  a  matter  of  convenience.  To  determine  the  sign  of  ys , 
however,  d23  must  be  calculated  using  a  positive  and  a  negative  ys  and  compared  to  the 
actual  distance  between  stations  two  and  three.  That  is, 

d23=^|{x2-x^f +iy2-y2)  and  ^^23  =  V(^2 -^3)^ +  (^2  “^3 )  •  (3-19) 

CLIMAREF  uses  the  ys  value  yielding  the  d2s  that  matches  the  actual  distance  between 
stations  two  and  three.  Once  again,  NS  and  NIK  values  for  the  three  stations  are 
normalized  to  sea  level  and  used  as  the  z  coordinates. 

The  equation  describing  the  plane  containing  points  Sin,  S2n,  and  Ssn  is 

A(x-Xi)-l-B(y-y,)-l-C(z-Zi)  =0  (3-20) 

where  A,  B,  and  C  are  the  direction  numbers  of  a  line  perpendicular  to  the  plane.  Two 
vectors  lying  in  the  surface  of  the  plane  are: 
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“12  =  (^2  -  )'  +  (^2  -  Jl ) i  +  (Zi  -  Zi  )k 


(3-21) 


and  Mi3  =  (xj  -  )i  +  ();3  -  )j  +  (23  -  Zi  )k 

and  a  vector  perpendicular  to  the  plane  is  given  by  m,2  ®  M13 ,  so  that 

A  =  (>2  -  >'1  )(Z3  -  Zi )  -  (Z2  -  )(>'3  -  3^1 ) 

B  =  {Z2-Zi){x^-Xi)-{X2-Xi)iZ3-Zi) 

C  =  (X2-  X,  )(y3  -  y, )  -  (^2  -  Ji  )(^3  -  ^1 )  • 

Setting 

D  =  -Axj  -  Byj  -  Czi  , 

the  plane  is  defined  by 

Ax  By  -l-  Cz  D  —  0  . 

Since  the  interpolated  N-value  is  the  point  on  the  plane  directly  above  point  r, 

^  interp  • 


(3-22) 


(3-23) 


(3-24) 


(3-25) 


(3-26) 


As  before,  CLIMAREF  de-normalizes  the  interpolated  values  to  the  interpolated  radar 
elevation  to  complete  the  interpolation  process. 

Once  the  refractivity  parameters  are  in  hand,  the  profile  is  constructed.  To  do  this, 
the  atmosphere  is  sliced  horizontally  into  layers  dh  thick.  Then  the  refractivity  N(h)  is 
calculated  for  every  layer  boundary  from  the  surface  up  to  a  maximum  altitude  five 
kilometers  higher  than  the  radar  height  or  the  target  height,  whichever  is  higher.  The  five 
kilometers  is  a  buffer  to  provide  for  circumstances  in  which  the  trace  itself  actually  goes 
higher  than  either  endpoint  (This  could  occur  in  a  duct). 

If  the  user  selects  not  to  include  ducting  in  the  profile,  the  calculations  are  simple, 
again  using  equations  2-8  and  2-9: 
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N(h)  =  Ns  exp{-c^ (h-hs)}. 


(2-8) 


where 


c,  =  In 


N 


s  _ 


=  ln- 


Ns+AN 


(2-9) 


If  ducting  is  to  be  included,  however,  every  unique  level  of  the  duct  must  be  taken  into 


account  when  building  the  profile. 


(a) 


(b) 


Surface-Based  Duct  Parameters 

GMSB 

M-unit  gradient 

OHSB 

optimum  height 

MTSB 

thickness 

MDSB 

M-unit  deficit 

MFSB 

trapping  frequency 

PSB12 

1 200Z  %  occur 

PSBOO 

OOOOZ  %  occur 

Elevated  Duct  Parameters 

GMEL 

M-unit  gradient 

OHEL 

optimum  height 

MTEL 

thickness 

MDEL 

M-unit  deficit 

MFEL 

trapping  frequency 

PEL12 

1 200Z  %  occur 

PELOO 

OOOOZ  %  occur 

P2EL 

possibility  of  >1 
elevated  duct 

PSBEL 

probability  of 
surface-based 

and  elevated  ducts 

Figure  3.9  Duct  parameters:  (a)  Surface  duct,  (b)  Elevated  duct  (figure  adapted  from 
Patterson,  1987) 


Figure  3.9  illustrates  how  each  duct  is  built  up  from  the  parameters  in  the  database. 


Before  the  N(h)  values  are  computed,  CLIMAREF  calculates  the  profile  height,  ph,  and 
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refractivity,  pM,  at  each  of  the  significant  levels  in  the  duct  using  the  following  formulas 
derived  from  the  figure: 

Surface  Duct: 


surface: 

phi  =0  pMi=:  NS 

(3-27) 

opt  coupling  ht: 

ph^  =  OHSB  pM^  =  NS  +  GMSB  ■  OHSB 

(3-28) 

top  of  duct: 

ph^  =  MTSB  pM^  =  pM^  -  MDSB 

(3-29) 

1  km  above  duct: 

ph^  =  ph^  +  \km  pM^  =  pM^  +  (MIK  -  NS) 

(3-30) 

Elevated  Duct: 

surface: 

phi  =  0  pMi  =  NS 

(3-31) 

opt  coupling  ht: 

ph2  =  OHEL  pM  2  =  NS  +  GMEL  ■  OREL 

(3-32) 

top  of  duct: 

MDEL 

phj  -  OHEL  -  ^ +  MTEL  pM.  -  pM,  -  MDEL 

GMEL 

(3-33) 

1  km  above  duct: 

ph^  =  plh,  +  \km  pM^  =  pMj  +  {MIK  -  NS) 

(3-34) 

All  the  heights  given  here  are  measured  with  respect  the  surface. 

Once  the  program  arrives  at  the  M-values,  it  converts  them  to  N- values  using 

N{h)=  M{h)-\5eh,  (3-35) 

where  h  is  the  height  above  the  surface,  following  the  convention  of  the  database 
(Patterson,  1987,  14)  (see  also  eqs.  2-14,  2-15) .  Above  the  top  of  the  duct,  CLIMAREF 
calculates  the  profde  using  pNs,  pN4  as  NS  and  NIK ,  respectively,  in  equations  2-8  and 
2-9.  With  the  profile  set  up,  the  stage  is  set  for  raytracing. 

Raytracing.  Raytracing,  as  the  term  is  used  generally  in  this  paper,  is  the  process 
of  computing  the  single  most  direct  (or  shortest)  path  from  the  radar  to  the  target  as  it  is 
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constrained  by  the  prevailing  refractive  effects  of  the  atmosphere.  This  path,  due  to 
refraction,  will  not  be  a  straight  line;  neither  will  it  be  a  function  with  a  closed  solution. 
Hence,  the  only  way  this  path  can  be  computed,  in  the  general  case,  is  by  breaking  the 
atmosphere  into  exceedingly  thin  concentric  shells  around  the  earth,  determining  the 
index  of  refraction  for  each  layer,  launching  a  ray  into  space  at  some  initial  takeoff  angle 
from  the  radar  antenna,  and  calculating  how  much  the  ray  bends  as  it  passes  through  each 
successive  layer.  The  fulcrum  upon  which  the  whole  process  turns  is  Snell’s  Law,  as 
shall  be  seen. 

Unfortunately,  the  only  way  of  putting  the  endpoint  of  the  ray  squarely  upon  the 
target  is  the  process  of  estimation,  approximation  and  iteration  described  in  the  first 
paragraphs  of  this  chapter.  The  details  will  be  described  presently.  First,  however,  it  is 
important  to  define  the  fundamental  types  of  paths  that  can  be  produced  by  such  a 
process. 

When  no  ducting  is  present,  there  are  three  path  types.  The  simplest  occurs  when 
the  radar  elevation  angle  is  positive  (Figure  3.10a).  The  energy  travels  upwards  with  a 
gradual  bend  toward  the  earth  until  it  reaches  the  target.  On  the  other  hand,  if  a  ray  is 
launched  at  a  negative  angle  it  will  travel  towards  the  earth,  bending  gradually  towards 
the  earth  and  eventually  either  strike  the  surface  (and  be  reflected  or  absorbed)  or  pass 
close  to  the  surface  at  some  minimum  height  (the  tangent  point)  only  to  continue  on  into 
space.  From  this  scenario,  two  path  types  arise:  the  short  path  (Figure  3. 10b),  on  which 
the  ray  reaches  the  target  before  hitting  the  earth  or  passing  through  the  tangent  point;  and 
the  long  path  (Figure  3.10c),  on  which  the  ray  only  reaches  the  target  after  it  has  passed 
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through  the  tangent  point.  The  dueted  path,  will  either  be  a  distorted  variant  of  one  of  the 
three  non-ducted  paths,  or  it  will  be  a  path  fully  contained  within  the  duct.  This  last  path 
type  will  only  occur  if  the  radar  and  the  target  are  in  the  duct  together. 


Figure  3.10  Basic  path  types:  (a)  upward,  (b)  downward-short,  (c)  downward-long 
(figure  adapted  from  Abel,  1982:  5) 


Of  course,  these  paths  are  only  realized  if  CLBSdAREF  can  put  the  ray  endpoint  on 
the  target.  Mathematically  the  problem  is  this:  For  a  given  radar  height,  target  height, 
target  range  and  refractive  profile,  a  ray  launched  at  an  initial  elevation  (takeoff)  angle, 

Oo,  will  cause  the  ray  to  pass  through  the  target  height  at  an  estimated  subtense,  I3e((Xo) 
(Figure  3.1).  That  subtense  will  be  correct  when, 

(3,(ao)  =  p,  ,  or  p,(ao)-|3,  =/(ao)  =  0  (3-36) 
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where  |3t  is  the  target  subtense  calculated  from  the  target  range,  target  height,  earth’s 
radius,  and  the  surface  elevation.  The  algorithm  used  for  this  thesis  project  is  taken  from 
USAFETAC  Technical  Report,  TN-82/005,  “The  Theory  and  Use  of  a  Raytracing  Model 
Developed  at  USAFETAC,”  by  Capt  Michael  D.  Abel,  et.  al.,  published  in  1982.  Abel’s 
method  is  an  example  of  a  simple  iterative  technique  for  approximating  a  root  known  as 
the  Newton-Raphson  method  of  numerical  approximation.  The  Newton  method  is 
summarized: 


(3-37) 


In  the  simplest  cases,/’,  is  known.  However,  in  this  case,  because/fcxo)  cannot  be  put  in 
closed  form,  it  must  be  numerically  approximated  on  the  fly  as  will  be  shown. 

Using  the  Newton-Raphson  Method,  then,  the  first  step  is  to  determine  an  initial 
estimate  for  the  takeoff  angle.  To  do  this  CLIMAREF  uses  4/3-earth  standard  refraction, 
an  idea  suggested  by  Abel  (1982:  21)  though  unaccompanied  by  any  algorithms.  Hence, 
the  algorithm  for  making  this  initial  guess  is  derived  here.  Enlarging  the  radius  of  the 
earth  by  1/3  (Figure  3.1 1)  and  keeping  the  radar  height,  target  height  (for  this  derivation 
surface  elevation  is  disregarded  for  simplicity),  and  surface  range  constant  under  4/3- 
earth  refractive  conditions  causes  the  radio  path  to  flatten  into  a  line,  and  the  radar-to- 
target  subtense  to  decrease  by  1/4,  i.e. 

CL  3 

P4/3  =  Pr^=PrT-  (3-38) 


The  4/3 -earth  model  and  the  Law  of  Cosines  show  the  apparent  range  (path  length)  to  be 
=  7 (fl4/3  +  RElf  +  (a4/3  -1-  THTf  -  +  REL){a^,^  +  THT)  cos  P4/3  ,  (3-39) 
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allowing  the  law  of  sines  to  be  used  to  find  our  initial  estimate: 


ttn 


cos 


_if  (<^4/3  +  THT)  sin  P4/3 


TRG 


APP 


(3-40) 


Figure  3.1 1  4/3  earth  estimation  of  initial  takeoff  angle;  (a)  earth  actual  size,  refracted  ray 
path  (b)  earth  enlarged  to  4/3  actual,  ray  path  subsequently  flattened 

CLIMAREF  now  uses  this  initial  estimate  of  the  ray  takeoff  angle  to  begin  the 
first  raytrace.  Since  the  basic  raytrace  process  is  discussed  in  some  detail  in  chapter  two, 
only  the  equations  are  given  here  using  the  variable  names  used  in  CLIMAREF. 
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Following  those  we  shall  examine  the  somewhat  complex  estimation  and  approximation 
techniques. 

Beginning  at  the  radar  and  ascending  or  descending  to  the  target,  the  bending  of 
the  ray  must  be  calculated  over  every  layer.  As  the  ray  approaches  a  layer,  the  known 
parameters  are  the  initial  elevation,  tti,  the  thickness  of  the  layer,  dh,  and  the 
refractivities,  Ni  and  N2 ,  at  boundary  1  and  boundary  2  of  the  layer.  The  program  will 
compute  the  elevation  angle  of  the  ray  as  it  approaches  boundary  2,  CL2,  the  amount  of 
bending,  'Fi,  through  the  layer,  and  finally,  the  subtense,  py,  of  the  ray  path  through  the 


layer.  These  are  calculated  like  so  (eq  3-41,  3-42,  3-43  are  from  Abel,  1982, 15-16): 


dir-dh\ 


cosa2=  1-(-(Ai-A^2)x10  - - T 

L  ^0 


(3-41) 


(where  dir  is  the  direction  of  the  ray  through  the  layers:  1  indicates  up  -1  indicates  down) 


Figure  3.12  Basic  raytracing  geometry 
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2(iVi-iV2)xlO^ 
tana,  +tana9 


(3-42) 


Pi  —  'Fj  +  ctj  tti 


(3-43) 


TRGappi  is  calculated  using  the  Law  of  Cosines  as, 


TRGappi  =  -^(^0  +/i,)^  +  (ro  +  /!2)^  -2(^0  -l-/i,)(ro  4-/12)  cos  Pi 


(3-44) 


Note  that  this  length  is  approximated  as  the  secant  of  the  actual  path  over  the  layer.  Next, 
Pi  is  accumulated  into  the  total  estimated  subtense,  Pe,  and  TRGappi  is  accumulated  into 
the  total  apparent  range,  TRGapp-  The  next  step  is  to  calculate  5pi/5ao  over  the  layer: 


5Pi  tantto  tantto  'Fitantto  [1  1 

dtto  tana,  tana2  tanaj -i-tana2  i^cosaj  sina,  cosa2sina 


(3-45) 


This  result  is  accumulated  and  used  to  estimate  the  next  initial  takeoff  angle.  Finally, 
CLIMAREF  updates  the  heights,  boundary  1  elevation  angle,  and  the  refract! vity  values, 
then  loops  back  and  performs  the  calculations  for  the  next  layer. 

In  the  case  of  a  downward  long  path  or  a  ducted  path  the  trace  will  pass  through  a 
height  extrema,  i.e.  a  point  at  which  the  ray  is  tangent  to  the  earth’s  surface.  When  cosa2 
(eq.  3-41)  comes  out  to  be  greater  than  one,  CLIMAREF  knows  an  extrema  has  been  hit 
and  that  special  treatment  is  required  (Abel,  1982:  23-29).  First,  the  extrema  height, 
which  generally  will  not  fall  on  one  of  the  predetermined  height  levels,  must  be 
calculated  as  (Abel,  1982:  23-24): 


\b\-4b^-ac 

h,n  =  hy+  dir  ■  Az  —  1\+  dir  ■ - 

A 


(3-46) 
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where, 


1  5n 

rg  +  /zj  8h 

1  Sn'^ 
2[rQ  +  hi  Sh  j 


(3-47) 

(3-48) 


C  cosaj  .  (3-49) 

Also,  6Pi/5oco  is  calculated  differently  for  the  tracing  steps  on  either  side  of  the  extrema. 
The  value  prior  to  the  extrema  is  (Abel,  1982, 29): 

50) 

The  value  for  the  step  following  the  extrema,  by  symmetry,  must  be  the  same,  although 
Abel  has  two  separate  formulas  for  the  two  steps.  (Note:  The  first  term  in  equation  44  in 
Abel  should  have  a  negative  sign  in  front  of  it,  and  the  second  formula  does  not  result  in 
the  same  value  as  the  first  as  it  should  for  reasons  undiscovered  by  this  author.)  Once  the 
extrema  height  has  been  established,  the  rest  of  the  layer  calculations  proceed  as  usual. 

Generally  speaking,  the  trace  should  stop  when  it  reaches  the  target  height,  and  for 
the  upward  path  this  simple  logic  is  all  that  is  required.  However,  for  downward  paths 
and  paths  where  ducting  is  present,  a  more  involved  set  of  rules  is  required  to  tell  the 
program  when  to  stop  a  given  trace  and  reestimate  if  necessary.  To  wit,  nine  unique 
stopping  cases  have  been  defined:  The  trace  will  stop  if: 

Case  1 :  Radar  and  target  not  in  duct  together,  Oo  positive,  trace  is  ducting 
Case  2:  Radar  and  target  not  in  duct  together,  Oq  positive,  target  height  reached 


-tantto  2no?bsintto 

- + - 7 - - 

tana, 

tana,  —  +  rr.+h) 

\y 


'F,  tanaQ 
sin^  a. 
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OR 


Radar  and  target  are  in  duct  together,  cxo  positive,  target  height  reached,  then 
trace  leaves  duct 

Case  3:  Radar  and  target  not  in  duct  together,  Oo  negative,  trace  is  ducting 
Case  4:  Radar  and  target  not  in  duct  together,  (Xo  negative,  target  is  higher  than  radar, 
target  height  reached 

OR 

Radar  and  target  are  in  duct  together,  ocq  negative,  target  is  higher  than  radar, 
trace  leaves  duct,  target  height  reached 

Case  5:  Radar  and  target  not  in  duct  together,  Oo  negative,  target  is  lower  than  radar, 
target  height  reached  twice  (long  path) 

OR 

Radar  and  target  not  in  duct  together,  ckq  negative,  target  is  lower  than  radar, 
target  height  reached  once  and  trace  reached  the  surface  of  the  earth  (short  path) 

OR 

Radar  and  target  are  in  duct  together,  oto  negative,  target  is  lower  than  radar, 
trace  leaves  duct,  target  height  reached  twice  (long  path) 

OR 

Radar  and  target  are  in  duct  together,  oto  negative,  target  is  lower  than  radar, 
trace  leaves  duct,  target  height  reached  once  and  trace  reached  the  surface  of  the 
earth  (short  path) 
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Case  6:  Radar  and  target  are  in  duct  together,  trace  remains  in  duct,  trace  endpoint  is 
beyond  target,  target  height  reached 

OR 

Radar  and  target  are  in  duct  together,  oco  positive,  trace  passes  through  target 
eight,  reaches  maxima,  trace  reaches  surface  of  earth  on  downslope 
Case  7:  Oo  negative,  trace  reached  surface  of  the  earth,  and  no  other  stop  situation  exists 
Case  8:  Radar  and  target  not  in  duct  together,  Oo  negative,  target  is  lower  than  radar, 
trace  reaches  height  minima  at  a  point  higher  than  target 
Case  9:  Radar  and  target  are  in  duct  together,  trace  is  ducting  without  reaching  target 
height  (i.e.  amplitude  of  oscillating  duct  path  not  large  enough  to  reach  target) 
Some  stop  cases  will  trigger  CLMAREF  to  check  for  endpoint  on  target,  some 
will  not.  If  the  endpoint  of  the  ray  is  not  on  target,  a  new  estimate  for  cxo  is  computed. 
How  that  is  done  depends  directly  on  the  reason  for  stopping  the  previous  trace. 

If  the  trace  is  stopped  because  of  cases  2, 4,  5  or  6,  the  program  checks  to  see  if 
the  ray  subtense  matches  the  target  subtense  for  any  of  the  points  where  the  trace  passed 
through  the  target  height.  Here  “points”  is  plural  because  for  both  the  downward  path 
and  the  ducted  path,  the  trace  may  pass  through  the  target  height  more  than  once.  The 
tolerance  for  a  subtense  match  is  set  to  ±  10'^  radians,  which  corresponds  to  a  range 
tolerance  of  about  ±  6  meters  at  the  earth’s  surface.  For  simple,  non-ducting  scenarios, 
CLIMAREF  usually  beats  this  tolerance  in  two  to  four  iterations.  The  maximum 
allowable  iterations  is  currently  twenty.  There  are  cases  in  which  the  program  needs  the 
extra  iterations  to  zero  in  on  a  radio  hole  or  overcome  distortion  in  the  trace  caused  by  a 
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duct.  For  some  scenarios,  ducting  causes  enough  ambiguity  to  completely  baffle  the 
estimating  routines  and  the  iteration  limit  is  reached  forcing  an  error  message.  The  next 
two  chapters  cover  these  limitations  in  more  detail. 

If  the  trace  misses  the  target,  oto  is  used  in  the  next  estimation  and  either 
discarded  or,  if  it  puts  the  trace  closer  to  the  target  than  any  previous  angle,  stored  as  a 
limiting  takeoff  angle  for  either  the  high  side  or  low  side  of  the  target  subtense.  This  puts 
a  bound  on  any  future  estimated  takeoff  angles  to  guard  against  wild  estimations,  which, 
under  certain  ducting  conditions,  can  easily  occur.  Unfortunately,  it  is  not  a  perfect 
defense  against  a  diverging  approximation. 

Cases  1,  3,  7,  8,  and  9  do  not  require  a  target-reached  check  since  they  represent 
traces  aborted  for  reasons  other  than  target-height-reached. 

If  CLIMAREF  determines  that  the  ray  endpoint  is  not  yet  on  the  target,  it 
estimates  a  new  Oo  based  on  everything  it  knows  about  the  situation.  Generally,  the 
Newton-Raphson  method  will  be  used  to  make  the  new  estimation,  that  is,  from  eq.  3-38 


(Abel,  1982:  27), 


a'o  =  tto  + 


5p  /  5ao 


(3-51) 


where  a’o  is  the  new  takeoff  angle  estimate,  Oo  is  the  old  estimate,  pt  is  the  target 
subtense,  Pe  is  the  subtense  of  the  endpoint  of  the  last  raytrace  and  6p/6cxo  is  the 
accumulated  estimate  of  how  pe  changes  with  Oo  as  described  above. 
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Depending  on  the  stopping  case,  this  estimate  may  be  augmented  or  replaced  by 
other  estimating  techniques  designed  to  keep  the  tracing  out  of  trouble.  Here  are  the 
estimation  procedures  for  each  stopping  case: 

Case  1:  a)  If  every  trace  so  far  has  been  in  the  duct,  multiply  last  takeoff  angle  by  1.5  to 
break  out  of  the  duct 

b)  If  a  previous  trace  was  not  ducted,  estimate  Oo  by  splitting  the  difference 
between  the  lowest  out-of-duct  takeoff  angle  and  the  highest  ducted  takeoff  angle 
in  order  to  get  out  of  the  duct  without  going  too  high.  Note:  The  out-of-duct 
trace  was  too  high  since  otherwise  the  estimator  would  not  have  gone  lower 
causing  the  ray  to  be  ducted  again  (Figure  3.13). 


target 


Figure  3.13  Takeoff  angle  estimation  (flat  earth  representation)  -  Case  1 

Case  2:  Standard  Newton-Raphson  estimation 

Case  3:  a)  If  every  trace  so  far  has  been  in  the  duct,  multiply  last  takeoff  angle  by  1.5  to 
break  out  of  the  duct 
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b)  If  a  previous  trace  was  not  ducted,  estimate  Oo  by  splitting  the  difference 
between  the  highest  out-of-duct  takeoff  angle  and  the  lowest  ducted  takeoff  angle 
in  order  to  get  out  of  the  duct  without  going  too  low.  Note:  Similar  to  Case  1, 
except  Oo  is  negative. 

Case  4:  Standard  Newton-Raphson  estimation 
Case  5:  Standard  Newton-Raphson  estimation 
Case  6:  a)  Standard  Newton-Raphson  estimation 

b)  If  the  resulting  takeoff  angle  is  lower  than  the  last  angle  that  was  too  low  to 
allow  the  ducted  ray  to  reach  the  target  height,  multiply  the  estimated  angle  by 
1.2  to  increase  its  magnitude  (Figure  3.14). 


Figure  3.14  Takeoff  angle  estimation  -  Case  6 

Case  7:  a)  If  all  traces  have  hit  the  ground  and  target  is  higher  than  radar,  force  a 

positive  elevation  angle  0.25  the  magnitude  of  the  last  takeoff  angle.  This  will 
get  the  trace  out  of  the  dirt  and  force  the  approximation  to  converge  down  to  the 
proper  angle 
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b)  If  all  traces  have  hit  the  ground  and  the  target  is  lower  than  the  radar,  multiply 
the  last  Oo  by  0.75  to  eventually  get  it  out  of  the  dirt. 

c)  If  a  previous  trace  did  not  hit  the  ground,  split  the  difference  between  the 
highest  takeoff  angle  that  hit  the  ground  and  the  lowest  angle  that  did  not  hit  the 
ground.  The  logic  is  similar  to  Case  1,  b. 

Case  8:  a)  If  no  trace  has  yet  reached  down  to  the  target  height,  multiply  last  (Xo  by  1.2 
b)  If  a  trace  has  reached  down  to  the  target  height,  split  the  difference  between 
highest  trace  that  reached  target  height  and  the  lowest  that  did  not  (Figure  3.15). 


trace  does  reach  down  to  target  height 
Figure  3.15  Takeoff  angle  estimation  -  Case  8 


Case  9:  Multiply  takeoff  angle,  Oo,  by  1.2  to  allow  trace  to  reach  target  height 
Note  that  in  any  case,  if  the  target  height  has  been  reached  more  than  once,  the  (5e  used  in 
the  estimation  must  be  the  one  nearest  to  the  target,  naturally.  All  these  conditional 
adjustments  to  the  standard  approximation  are  necessary  to  cover  situations  that  are 
outside  the  simple  paradigm  assumed  by  the  Newton-Raphson  method. 
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Once  a  new  do  has  been  estimated,  CLIMAREF  compares  it  to  the  previous 
takeoff  angles  that  put  the  trace  endpoint  nearest  the  target  —  one  for  the  high  side  and 
one  for  the  low.  If  the  new  angle  is  outside  these  limits,  it  is  recalculated  as  the  bisection 
of  the  two  limiting  angles.  This  bisection  method  is  simpler  and  not  as  efficient  as 
Newton-Raphson,  but  because  ducting  effects  sometimes  cause  the  latter  method  to 
malfunction,  bisection  is  useful  as  a  backup  and  a  check.  When  the  new  takeoff  angle  is 
finalized,  a  new  trace  is  performed.  This  process  will  repeat  itself  until  the  trace  endpoint 
lands  squarely  on  the  target,  indicating  the  most  direct  path  to  the  target  has  been  found. 

When  the  ra5^racing  is  complete  with  the  traee  endpoint  on  the  target,  all  that  is 
left  to  do  is  calculate  the  range  and  height  error,  whieh  are  the  fundamental  output  of  the 
model.  First,  the  geometric  elevation  angle,  0%  (the  elevation  angle  at  the  radar  of  the 
straight  line  path  to  the  target),  of  the  target  is  computed. 


dg  =±asin 


(t,  +  THT)  sin  p A  n 


(3-52) 


To  determine  the  sign  of  ag  we  must  use  it  in  the  Law  of  Cosines  to  calculate  the  radius 
to  the  target,  then  compare  that  with  the  known  target  radius,  ro+THT. 

+  TRG^  - 2rJRGcos(ag  +n/ 2)  (3-53) 

Whichever  ttg  causes  r,  to  come  closest  to  r,=ro+THT,  is  the  one  used  for  the  remaining 
calculations. 

Next,  the  app^lrent  (or  measured)  target  height,  THTapp,  is  computed. 

THTj^PP  =  7?o  +  TRGIpp  -  IrJRG^pp  cos(ao  +  7i  /  2)  -  (3-54) 

With  these  parameters  in  hand,  CLIMAREF  can  calculate  the  errors: 
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(3-55) 


THTERR  =  THT^pp  -  THT 


THTERR 

THTERR,,,-  xlOO 

(3-56) 

TRGERR  =  TRGj,pp  -  TRG 

(3-57) 

TRGERR 

77?G£7?/?iaq  =  xlOO 

100 

(3-58) 

Actually,  the  absolute  error  is  probably  more  important  than  percent  error  for  air  traffic 
considerations,  but  both  are  calculated  as  a  matter  of  course. 

Display  Data.  Of  course  the  last  step  is  to  display  this  data  for  the  user.  Actual 
height  and  range,  apparent  height  and  range,  and  the  height  and  range  errors  are  alt 
displayed.  In  addition,  the  actual  path,  the  geometric  (straight-line)  path,  and  the  apparent 
path  are  plotted  in  one  of  two  representations:  curved  earth  (Figure  3.16)  or  flat  earth 
(Figure  3.17).  In  neither  representation  does  the  height  axis  match  the  range  axis.  If  the 
plots  were  to  scale,  no  distinction  could  be  made  between  the  various  curves.  They 
would  be  so  close  as  to  be  on  top  of  each  other. 
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So,  using  Snell’s  Law-based  raytracing  techniques  driven  by  climatologically 
predicted  refractive  profiles,  CLMAREF  predicts  radar  height  and  range  error  for  any 
time  of  year  at  nearly  any  location  in  the  world.  The  following  chapters  contain  the 
results  and  analysis  of  a  variety  of  proof  tests  demonstrating  the  usefulness  of 
CLIMAREF  for  its  intended  purpose.  Additionally,  some  of  the  program’s  capabilities 
and  limitations  are  presented. 
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IV.  Data  and  Analysis 


The  CLIMAREF  model  was  developed  for  this  thesis  as  a  new  approach  to 
height  error  prediction,  replacing  the  old  four-thirds  earth  and  exponential  standard 
atmosphere  techniques.  The  purpose  of  this  chapter  is  to  demonstrate  the  model’s 
validity,  to  present  sample  data  as  evidence  of  the  model’s  usefulness,  and  to  show  where 
the  algorithm  breaks  down.  First,  CLIMARFF  output  is  compared  to  a  height-range- 
angle  chart  derived  in  Blake  (1986:  185)  and  to  RAYS,  a  submodule  of  the  Navy’s 
FRFPS  software.  Next,  a  validation  test  of  the  interpolation  routines  demonstrates  the 
accuracy  of  the  geographical  calculations  involved.  And  then,  the  problem  of  a  minimum 
initial  ray  elevation  angle  is  explored.  After  that,  a  survey  of  climate  variation  in  terms  of 
height  error  is  presented.  And,  finally,  an  inside  look  at  CLIMARFF’ s  workings  provides 
insight  into  tracing  scenarios  where  the  estimation  algorithm  breaks  down.  Chapter  5 
contains  conclusions  suggested  by  the  results  presented  in  this  chapter. 

Validation:  CLIMAREF  vs.  Blake  and  CLIMAREF  vs.  RAYS 

The  best  way  to  validate  CLIMAREF  would  be  to  compare  its  results  with 
statistics  derived  from  an  actual  radar-target-atmosphere  system.  Such  a  comparison  was 
impossible  in  this  case  due  to  time  and  resource  constraints.  As  a  substitute,  the 
performance  of  the  raytracing  algorithms  used  in  CLIMAREF  was  compared  to  the 
height-range-angle  chart  in  Blake  (1986:  185)  and  to  FRFPS  RAYS. 
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The  height-range-angle  chart  in  Blake  allows  determination  of  the  proper  takeoff 
angle  for  a  given  target  height  and  range.  It  is  based  on  a  closed-form  geometrical  optics 
solution.  The  height  selected  was  30,000  feet,  and  the  range,  80  nautical  miles.  The 
required  takeoff  angle  given  by  Blake  was  3.0  degrees  as  close  as  it  could  be  read  on  the 
chart.  CLMAREF  was  modified  to  display  the  final  takeoff  angle  and  run  using  the 
standard  atmosphere  option.  The  final  takeoff  angle  given  by  CLIMAREF  was  3.0027 
degrees. 

The  second  test  was  also  fairly  simple. .  RAYS  “traces  the  paths,  in  height  and 
range,  of  electromagnetic  rays  based  upon  a  linearly  segmented  refractivity-versus- 
altitude  profile(s)...”  The  software  was  developed  and  is  used  by  the  Naval  Command 
Control  and  Ocean  Surveillance  Center.  (Patterson,  et.  al,  1994:  35).  To  make 
comparisons  with  RAYS,  a  program  called  RAYSD  was  developed  using  the  same 
tracing  mathematics  used  in  CLIMAREF.  RAYSD  differs  from  CLIMAREF  in  that  no 
target  is  involved.  The  user  selects  a  radar  height  and  a  ray  takeoff  angle  and  the  ray  is 
traced  to  a  maximum  height  or  range,  whichever  comes  first.  The  program  simply 
computes  the  path  of  a  given  ray  through  a  given  atmosphere.  EREPS  RAYS  does  much 
the  same  thing.  Unfortunately,  RAYS  is  a  DOS  program  and  the  output  is  to  the  screen 
only  -  there  is  no  data  output.  Nevertheless,  using  a  screen  capture  program  and 
overlaying  the  RAYS  and  RAYSD  plots,  the  degree  of  correspondence  between  the  two 
models  can  be  reasonably  ascertained. 

The  test  was  performed  first  using  a  surface  duct  and  then  an  elevated  duct 
scenario.  The  elevation  and  M- values  in  Table  4.1  were  used  in  both  programs  to 
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construct  the  refractivity  profiles.  For  the  surface  duct  test,  the  radar  height  was  3  feet 
and  traces  were  launched  using  takeoff  angles:  0.05, 0.1 125, 0.175, 0.2375,  and  0.3 
degrees.  For  the  elevated  duct  test,  radar  height  was  2500  feet  and  the  ray  takeoff  angles 
were  -0.5,  -0.3,  -0.1,  0.1  and  0.3  degrees.  The  results  are  plotted  in  Figure  4.1  and  Figure 
4.2. 

In  both  scenarios,  the  RAYSD  plots  match  the  RAYS  plots  perfectly  as  well  as 
can  be  detected  given  the  plot  resolution.  There  is  no  discernible  difference  between  the 
curves. 


Table  4. 1  Height  and  Modified  lOR  Parameters  for  RAYSD-RAYS  Validation  Tests  — 
Surface  Duct,  and  Elevated  Duct  Respectively 


Altitude  Above  Surface 

Modified 

IOR,M 

Altitude  Above  Surface 

Modified 

IOR,M 

feet 

meters 

feet 

meters 

0 

0 

323 

0 

0 

330 

174 

53 

329 

2369 

722 

414 

463 

141 

320 

2740 

835 

406 

3743 

1141 

435 

1835 

499 
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Height  Above  Surface  (Feet) 


Figure  4.1  RAYSD  vs.  RAYS  discrepancies  —  Surface  Ducts 


Figure  4.2  RAYSD  vs.  RAYS  discrepancies  —  Elevated  Duct 
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Interpolation  Check 

Two  interpolation  routines,  referred  to  as  the  linear  algorithm  and  the  planar 
algorithm,  were  developed  as  a  part  of  CLIMAREF  to  find  the  refractivity  values  near  or 
at  the  radar  site  using  data  from  two  or  three  nearby  radiosonde  stations.  To  verify  the 
accuracy  of  these  routines,  a  test  was  performed. 

Linear  algorithm.  To  test  the  linear  routine,  two  radiosonde  stations,  Pittsburgh 


R 


Figure  4.3  Geometry  for  test  of  linear  interpolation 

and  Buffalo,  were  chosen,  with  a  radar  located  off  center  somewhere  in  between  (40°N, 
80° W)  so  that  the  three  sites  were  configured  as  shown  in  Figure  4.3.  The  surface 
distances,  in  kilometers,  between  the  stations  were  obtained  from  the  distance-finding 
algorithm  in  the  CLIMAREF  module,  CR4calcdist  (which  was  also  checked  for 
reasonable  accuracy  against  a  map).  Distances  p,  b,  and  h  were  found  using  the 
Pythagorean  Theorem:  p  =  155.42,  b  =  139.06,  and  h  =  52.79.  Since  NS  and  NIK  at 
Buffalo  are  314.2348  and  278.7399,  respectively,  and  NS  and  NIK  at  Pittsburgh  are 
313.9588  and  279.0745  respectively,  a  manual  interpolation  along  line  PB  at  point  O, 
yields: 
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139.06  31 3.9588 

294.48  313.9588-313.2348  ^ 

which  is  used  as  NS  at  the  radar  site.  Using  the  same  site  locations  in  CRVinterp,  the 
output  was  NS  =  314.05.  The  difference  is  negligible  for  our  purposes.  Similarly,  the 
manual  computation  yielded  NIK  =  278.92,  which  exactly  matches  the  CRVinterp  NIK 
value. 

Planar  algorithm.  To  test  the  planar  interpolation,  a  third  station,  Flint-Bishop, 
Michigan,  was  selected.  The  radar  and  the  other  two  stations  were  left  alone.  The 
distances  between  the  stations  and  the  radar  were  again  derived  from  CR4calcdist  and  are 
given  in  Figure  4.4. 


P 


Figure  4.4  Geometry  for  test  of  planar  interpolation 


This  interpolation  is  a  bit  more  complicated:  First,  ([),  6,  and  \|/  were  found  using  the  Law 
of  Cosines.  Next,  a  and  P  were  derived  from  ([),  0,  and  \\f.  Then,  c//=313.14,  ^2=94.86, 
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and  <i5=107.34  were  found  using  the  Law  of  Sines.  Using  these  values  and  knowing  from 
CLIMAREF  that  NS  at  Buffalo  and  Flint  are  313.9588  and  312.6295  respectively,  NS  at 
point  O  was  linearly  interpolated  as  313.65.  Using  that  value  and  the  value  of  at 
Pittsburgh,  314.2348,  the  NS  at  the  radar  site  was  interpolated  to  be  313.88.  CLIMAREF 
computed  the  same  value.  Similarly,  the  manual  method  and  CLIMAREF  yielded  the 
same  value  for  NIK  at  the  radar  site,  278.7437. 

These  two  tests,  then,  demonstrate  the  validity  of  the  interpolation  algorithms 
developed  for  CLIMAREF. 


Run-Time  Test 

A  test  was  performed  to  determine  the  amount  of  processing  time  CLIMAREF 
requires  in  its  non-compiled  MATLAB®  coded  form.  The  test  was  run  on  three  different 
computers,  a  66  MHz  486  IBM  PC  compatible,  a  166  MHz  Pentium,  and  a  SUN 
Microsystems  Sparcstation  20.  The  MATLAB®  tic  and  toe  commands  were  placed  in 
CRlmain  before  and  after,  respectively,  the  call  to  CR9raytrace  in  order  to  time  the 
tracing  only.  The  results  are  listed  in  Table  4.  2. 


Table  4.  2  Run-time  results  for  three  targets  in  standard  atmosphere  refraction 


Radar  El  (ft) 

Targ.  Ht  (ft) 

Target  Rng 

(Nmi) 

#  Iterations 

Req’d 

486  66 

Pentium  166 

Sparc  20 

30 

2000 

5 

3 

35.2 

2.5 

8.0 

30 

10,000 

30 

3 

246.2 

21.0 

66.0 

30 

40,000 

150 

3 

out  of  mem 

209.7 

456.0 
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Minimum  Elevation  Angle 


On  rare  occasions,  the  initial  elevation  angle  for  a  CLIMAREF  trace  is  so  small 
that  problems  occur.  If  the  angle  is  too  small  and  negative,  the  first  subtense  calculated 
blows  up  to  an  unrealistically  large  value.  If  it  is  too  small  and  positive,  either  the 
subtense  blows  up  or  the  estimation  routine  has  hard  time  finding  the  correct  angle. 
Although  CLIMAREF  handles  these  errors  gracefully  (if  not  completely  satisfactorily)  it 
is  important  to  understand  the  nature  of  the  problem.  Hence  a  test  was  performed. 

Small  Negative  Angles.  To  examine  the  raytracing  response  to  a  small  negative 
angle,  CLIMAREF  was  run  using  standard  atmosphere  refraction  and  various  target 
ranges  and  heights  until  the  small,  negative  angle  response  was  observed.  The  particular 
scenario  was  radar  elevation  =  0.3  km  (984  ft),  target  height  =  0.305  km  (1001  feet),  and 
target  range  =  10  km  (5.4  Nmi). 

Everything  went  wrong  on  the  first  tracing  step.  CLIMAREF  chose  an  initial 
elevation  angle,  Oq  =  -5.65350x10"^  radians.  That  is,  it  expected  to  trace  a  downward 
path  that  would  reach  a  tangent  point,  and  find  the  target  as  it  moved  away  from  the  earth. 
The  first  level  height,  hj,  was  of  course  the  radar  height.  The  second  level  height,  then 


was  /i2=0.299  km.  Since  cos(ai)=0.999999998401,  and  using 


f-  dir  ■  dh 

cosa2=  l  +  (iVi-A^2)xlO  - - 7- 

L  - 


COSttj 


(3-41) 


together  with  the  appropriate  values  for  Nj,  N2,  dir,  dh,  and  ro,  cos(a2)  was  calculated 
as  1 .00000011240.  To  the  raytracing  routine,  a  result  greater  than  one  indicates  an 


extrema  has  been  reached,  in  this  case  a  minima.  So  tta  is  set  equal  to  zero  and  the 


bending,  v|/i ,  over  the  layer  is  computed  using: 


Vi  = 


2(ni-«2)xl0^ 
tanttj  +tana2 


(2-19) 


Since  ai  was  so  small  and  a2  was  zero,  \|/i  blew  up.  In  turn.  Pi, 

Pi  =  \l/i  -I-  a2  -  ttj  ,  (2-20) 

blew  up  causing  the  extreme  results  described  above. 

The  problems  disappeared  when  the  target  was  raised  to  0.306  km  and  when  it 
was  lowered  to  0.300. 

Small  Positive  Angles.  When  the  target  was  at  0.306  km,  a  single  raytrace 
worked  fine  using  a  positive  angle  of  2.655  x  10'^.  Further  testing  was  precluded  for  lack 
of  time. 


Climate  variation 

The  primary  concern  of  the  project  sponsor,  NAIC,  is  finding  a  technique  for  ray 
path  prediction  and  height  error  prediction  that  is  better  than  the  traditional  four-thirds 
earth  or  standard  atmosphere  solutions.  To  demonstrate  the  added  value  of  CLIMAREF, 
a  test  was  conducted  in  which  the  ray  path  was  computed  for  a  sample  target  under 
varying  conditions.  The  height  error  was  calculated  for  each  scenario  and  plotted  with 
reference  to  the  standard  atmosphere  height  error,  which  was  804  feet. 

The  sample  target  and  radar  positions  were  held  constant  throughout  the  test. 

The  radar  was  30  feet  above  the  surface,  the  target  was  10,000  feet  above  the  surface,  and 
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the  two  were  60  nautical  miles  apart.  The  variables  were  month:  January,  April,  July,  and 
October;  duct  type:  no  ducting,  surface,  duct,  elevated  duct;  and  location:  seventeen 
locations  representing  a  variety  of  North  American  and  world  climate  tj^es  (see  tables 
below).  Fifteen  of  the  seventeen  locations  were  chosen  to  correspond  to  locations  used 
by  Bean  and  Dutton  in  Radio  Meteorology  (1966:  109-131).  Where  the  database  was 
lacking,  substitutions  were  made,  i.e.  Cape  Kennedy  was  used  instead  of  Cocoa,  FL; 
Howard  AB  was  used  instead  of  Balboa,  Panama;  Port  Hardy,  Canada  was  used  instead 
of  Tatoosh  Is,  WA;  and  Wheelus,  Libya  was  used  instead  of  Tripoli,  Libya.  All 
substitutes  are  close  in  climate  and  proximity  to  the  originals.  The  other  two  of  the 
seventeen  stations,  Jodhpur,  India  and  Fort  Lamy  in  the  Sahara  Desert,  were  chosen 
because  Bean  and  Dutton  noted  those  areas  as  having  the  largest  ranges  of  surface  N-units 
(1966:  102).  The  results  are  tabulated  in.  Table  4.3,  Table  4.4  and  Table  4.5  and  shown 
plotted  in  Figure  4.5  through  Figure  4.1 1.  Table  4.6  holds  some  statistics  derived  from 


the  climate  variation  test. 


Table  4.3  Height  Error  (feet)  —  No  Ducting 
{REL='iO  feet,  T//r= 10000  feet,  77?G=60  Nmi) 


Location 

JAN 

APR 

JUL 

OCT 

Portland,  ME 

677.3 

695.7 

989.1 

806.2 

Washington,  DC 

695.9 

990.2 

788.1 

Hatteras,  NC 

1257.6 

1006.9 

Cape  Kennedy,  FL 

917.9 

1224.3 

1010.0 

Miami,  FL 

882.8 

901.0 

1154.0 

1028.6 

Howard  AB,  Panama 

1064.4 

1047.5 

1084.3 

1084.0 

Brownsville,  TX 

826.3 

1133.3 

1222.5 

955.6 

Columbia,  MO 

695.1 

695.6 

917.8 

732.2 

Bismarck,  ND 

748.7 

638.7 

787.2 

675.6 

Denver,  CO 

634.4 

579.6 

745.9 

San  Diego,  CA 

931.3 

1019.7 

1463.0 

Oakland,  CA 

914.2 

1394.8 

1142.8 

Port  Hardy,  Canada 

769.7 

844.0 

806.9 

Churchill,  Canada 

878.9 

768.9 

806.3 

732.5 

Jodhpur,  India 

675.7 

675.8 

920.0 

824.0 

Fort  Lamy,  Chad 

747.4 

711.2 

1116.1 

1110.7 

Wheelus,  Libya 

806.0 

1263.3 

1713.0 

1165.0 

Table  4.4  Height  Error  (feet)  —  Surface  Ducting 
(REL=30  feet,  77/7’= 10000  feet,  77?G=60  Nmi) 


Location 

JAN 

APR 

JUL 

OCT 

Portland,  ME 

1179.7 

1032.7 

Washington,  DC 

867.6 

922.4 

1187.0 

999.3 

Hatteras,  NC 

1015.9 

dbeh 

1595.0 

1313.1 

Cape  Kennedy,  FL 

1141.5 

illdef 

1504.4 

1195.6 

Miami,  FL 

1468.8 

1555.6 

1364.0 

Howard  AB,  Panama 

1253.0 

1177.9 

Brownsville,  TX 

1313.7 

1421.1 

1430.7 

1157.2 

Columbia,  MO 

798.7 

939.0 

1106.8 

1077.1 

Bismarck,  ND 

904.5 

898.6 

1101.3 

817.4 

Denver,  CO 

810.3 

770.8 

1018.7 

799.0 

San  Diego,  CA 

1118.3 

1205.4 

1905.6 

1583.8 

Oakland,  CA 

1048.0 

1172.8 

1615.2 

1366.2 

Port  Hardy,  Canada 

949.2 

989.7 

1003.7 

969.3 

Churchill,  Canada 

1057.3 

911.9 

1021.4 

945.7 

Jodhpur,  India 

1453.3 

1858.2 

1029.7 

Fort  Lamy,  Chad 

923.9 

1590.4 

1493.1 

Wheelus,  Libya 

1641.8 

2256.0 

1478.8 
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Table  4.5  Height  Error  (feet)  —  Elevated  Ducting 
{REL=30  feet,  THT=l(mO  feet,  TRG=60  Nmi) 


Location 

JAN 

APR 

JUL 

OCT 

Portland,  ME 

714.4 

693.4 

832.3 

1032.6 

Washington,  DC 

690.8 

728.1 

809.9 

Hatteras,  NC 

862.9 

945.2 

Cape  Kennedy,  FL 

931.1 

1020.1 

1235.4 

Miami,  FL 

920.4 

996.0 

1109.5 

950.0 

Howard  AB,  Panama 

961.2 

920.5 

1002.9 

868.9 

Brownsville,  TX 

951.7 

1186.6 

1233.0 

944.2 

Columbia,  MO 

776.0 

783.1 

899.1 

720.6 

Bismarck,  ND 

667.7 

792.2 

761.4 

Denver,  CO 

700.0 

678.0 

701.4 

San  Diego,  CA 

1028.9 

1094.7 

1596.2 

1287.5 

Oakland,  CA 

887.3 

1004.0 

1502.0 

1162.2 

Port  Hardy,  Canada 

801.6 

762.6 

870.6 

769.5 

Churchill,  Canada 

915.0 

809.6 

801.8 

Jodhpur,  India 

530.8 

258.7 

887.0 

Fort  Lamy,  Chad 

-1973.1 

460.55 

1073.2 

101.2 

Wheelus,  Libya 

788.5 

1323.7 

Table  4.6  Height  Error  Statistics  (REL=30  feet,  THT=10000  feet,  TRG=60  Nmi) 

Note:  Std  Atm  Ht  Err  =  804. 1  feet 


units=feet 

Mean 

Std  Dev 

Dev  from 
Std  Atm 

No  duct 

222.6 

249.3 

1713.0 

Surface  duct 

299.5 

483.7 

2256.0 

■IBBI 

Elevated  duct 

860.6 

431.5 

435.1 

1643.9 

101.2 

Overall 

987.0 

358.3 

402.3 

2256.0 

101.2 
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Minimum  Elevation  Angle 


On  rare  occasions,  the  initial  elevation  angle  for  a  CLIMAREF  trace  is  so  small 
that  problems  occur.  If  the  angle  is  too  small  and  negative,  the  first  subtense  calculated 
blows  up  to  an  unrealistically  large  value.  If  it  is  too  small  and  positive,  either  the 
subtense  blows  up  or  the  estimation  routine  has  hard  time  finding  the  correct  angle. 
Although  CLIMAREF  handles  these  errors  gracefully  (if  not  completely  satisfactorily)  it 
is  important  to  understand  the  nature  of  the  problem.  Hence  a  test  was  performed. 

Small  Negative  Angles.  To  examine  the  raytracing  response  to  a  small  negative 
angle,  CLIMAREF  was  run  using  standard  atmosphere  refraction  and  various  target 
ranges  and  heights  until  the  small,  negative  angle  response  was  observed.  The  particular 
scenario  was  radar  elevation  =  0.3  km  (984  ft),  target  height  =  0.305  km  (1001  feet),  and 
target  range  =  10  km  (5.4  Nmi). 

Everything  went  wrong  on  the  first  tracing  step.  CLIMAREF  chose  an  initial 
elevation  angle,  Oq  =  -5.65350x10’^  radians.  That  is,  it  expected  to  trace  a  downward 
path  that  would  reach  a  tangent  point,  and  find  the  target  as  it  moved  away  from  the  earth. 
The  first  level  height,  hi,  was  of  course  the  radar  height.  The  second  level  height,  then 


was  h2=0.299  km.  Since  cos(ai)=0.999999998401,  and  using 


£  dir  ■  dh 

cosa2  =  l-t-(A|-A/’2)xlO  - - —  cosaj 

L  ^0  ^  -I 


(3-41) 


together  with  the  appropriate  values  for  Nj,  N2,  dir,  dh,  and  ro,  cos(a2)  was  calculated 
as  1.0000001 1240.  To  the  raytracing  routine,  a  result  greater  than  one  indicates  an 
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Figure  4.7  Climate  variation  —  Elevated  Ducting,  Seasonal  Comparisons 
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Figure  4.8  Climate  variation  —  January,  Ducting  Comparisons 
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Figure  4.9  Climate  variation  —  April,  Ducting  Comparisons 
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Figure  4.10  Climate  variation  —  July,  Ducting  Comparisons 
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Capability  of  Model  to  Handle  All  Types  of  Refractive  Effects 

The  convenient  assumption  in  developing  a  radar  beam  bending  model  would  be 
that  the  radar  is  on  the  ground  and  the  target  is  in  the  air,  well  above  any  ducting  effects 
and  within  the  line-of-sight  range  of  the  radar.  Unfortunately,  such  is  not  the  case.  A 
general  purpose  model  must  be  able  to  respond  reliably  to  any  conceivable  situation.  In 
this  vein,  the  next  test  presents  the  model  with  radar-target-atmospherics  combinations 
representing  most,  if  not  all,  possible  combinations  of  radar,  duct  and  target,  to  determine 
where  the  weak  points  are,  if  any  exist. 

To  perform  this  test,  a  specially  modified  version  of  the  driving  module, 
CRlmain_tst,  was  used,  and  the  variable  plotit  in  CRraytrace9  was  set  to  1  allowing  all 
estimate-traces  to  be  plotted  along  with  the  final  trace.  In  the  plots  that  follow,  the  solid 
line,  if  any,  represents  the  final  raytrace  to  the  target  (represented  by  the  star),  and  the 
dashed  lines  represent  estimate-raytraces  leading  up  to  the  final.  These  estimate-traces 
were  numbered  in  order  by  CLIMAREF,  but  in  many  cases  their  proximity  caused  the 
numbers  to  be  obscured.  Dashed  horizontal  lines  represent  the  boundaries  of  a  duct. 
Further,  all  plots  are  of  the  flat-earth  variety.  Since  the  superrefracted  rays  have  a  radius 
of  curvature  greater  than  that  of  the  earth’s  surface  the  standard  upward  paths  seem  to 
curve  away  from  the  “flat”  earth. 

The  figures  are  organized  by  duct  type  and  then  labeled  with  path  type  and  any 
other  distinguishing  characteristics.  Path  type  is  upward  path,  short  downward  path,  long 
downward  path,  or  ducted  path.  There  are  also  variations  on  these  and  at  times,  because 
of  ducting  distortion,  the  path  type  definitions  become  blurred. 
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Figure  4.12  No  duct -- upward  path 


Figure  4.13  No  duct  -  short,  downward  path  (target  reached  before  ray  hits  surface) 
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Figure  4.14  No  duct  -  short  downward  path  (target  reached  before  tangent  point) 


Figure  4.15  No  duct  —  long,  downward  path  (target  lower  than  radar) 
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Above  Surface  (feet)  Height  Above  Surfai 


Ground  Range  (Nmi) 

Figure  4. 19  Surface  duct  —  upward  path  (target  in  radio  hole) 


Figure  4.20  Surface  duct  —  ducted  path 


Figure  4.21  Surface  duct  —  short,  downward  path  (radar  above  duct,  target  inside  duct) 
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Figure  4.22  Surface  duct  —  downward  path  (radar  abv  duct,  target  in  duct,  in  radio  hole) 


Figure  4.23  Elevated  duct  —  upward  path  (radar  below  duct,  target  above  duct  and  far) 
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Figure  4.27  Elevated  duct  —  short,  downward  path  (radar  in  duct,  target  below  duct) 
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Figure  4.28  Elevated  duct  —  long,  downward  path  (radar  in  duct,  target  below  duct) 


Figure  4.29  Elevated  duct  —  upward  path  (radar  above  duct,  target  above  duct) 
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Figure  4.30  Elevated  duct  —  long,  downward  path  (radar  above  duct,  target  above  duct) 
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Figure  4.31  Elevated  duct  —  long,  downward  path  (radar  above  duct,  target  above  duct  in 
radio  hole) 
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Ground  Range  (Nmi) 


Figure  4.32  Elevated  duct  —  short,  downward  path  (radar  above  duct,  target  below  duct) 


Figure  4.33  Elevated  duct  —  short,  downward  path  (radar  above  duct,  target  below  duct) 
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Figure  4.34  Elevated  duct  —  short,  downward  path  (radar  above  duct,  target  in  duct) 


Figure  4.35  Elevated  duct  —  long,  downward  path  (radar  above  duct,  target  in  duct) 


Figure  4.36  Elevated  duct  —  long,  downward  path  (radar  above  duct,  target  in  duct) 


Figure  4.37  Elevated  duct  —  short,  downward  path  (rdr  abv  dct,  tgt  in  dct,  in  radio  hole) 


CLIMAREF  correctly  converged  on  the  initial  takeoff  angle  (to  within  the  10'^ 
radian  tolerance)  in  20  of  the  26  representative  scenarios  presented  in  Figure  4.12  through 
Figure  4.37.  There  were  no  problems  at  all  when  no  duct  was  present.  When  the  target 
was  over  the  horizon,  the  program  responded  with  an  error  message  to  that  effect.  Even 
when  ducting  was  present,  the  program  responded  correctly  to  many  odd  situations. 
Particularly,  for  the  situations  illustrated  in  Figure  4.19  and  Figure  4.22  CLIMAREF 
correctly  reported  the  target  was  in  a  radio  hole.  Although  some  situations  took  more 
iterations  to  find  the  target  than  others,  in  the  simple  cases  three  usually  sufficed. 

Unfortunately,  six  of  the  scenarios  baffled  CLIMAREF  and  it  timed  out  after 
attempting  the  maximum  30  estimations.  Careful  examination  of  these  cases  (Figure 
4.20,  Figure  4.26,  Figure  4.30,  Figure  4.31,  Figure  4.36,  and  Figure  4.37)  reveal  two 
distinguishing  characteristics  of  traces  that  failed:  1)  both  the  radar  and  the  target  were  in 
the  duct  together  and  the  path  was  a  ducted  path,  or  2)  the  path  was  a  long  downward  path 
in  which  the  rays  crossed  (apparently  as  an  effect  of  an  elevated  duct). 

The  immediate  reasons  for  these  two  types  of  breakdown  are  readily  discernible: 

In  a  duct  (at  least  in  the  mathematical  kind  used  in  CLIMAREF),  there  are  many  paths 
(hence  takeoff  angles)  that  will  lead  to  a  given  target.  Because  of  how  the  rays  are 
refracted,  rays  starting  from  a  variety  of  takeoff  angles  will,  at  times,  cross  causing  more 
ambiguity.  A  similar  mechanism  causes  the  rays  in  the  downward  path  through  the 
elevated  duct  to  cross,  again  resulting  in  ambiguity. 

Without  further  discussion,  a  fundamental  problem  is  apparent.  One  of  the 
limitations  of  geometrical  optics  (see  chapter  2)  is  that  adjacent  rays  must  remain 
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approximately  parallel  to  each  other  within  a  wavelength.  In  other  words,  in  a  valid 
raytrace,  rays  can  not  cross.  This  is  an  important  observation  that  is  treated  further  in  the 
next  chapter. 

Numerical  analysis  provides  further  insight  into  the  problem.  The  estimation 
algorithm  used  by  CLIMAREF  depends  upon  accurate  knowledge  of  8p/5o(o,  that  is,  how 
much  the  subtense  of  the  ray  endpoint  changes  for  a  given  change  in  the  takeoff  angle. 
This  derivative  must  be  such  that  the  estimation  can  start  from  an  initial  guess  (four- 
thirds  earth  in  this  case)  and  find  the  takeoff  angle  that  will  result  in  the  correct  subtense. 
Figure  4.38  through  Figure  4.41  were  created  using  the  RAYSD  test  routine  to  plot  traces 
for  several  different  angles  under  various  circumstances.  In  each  figure,  one  or  more 
horizontal  lines  were  plotted  to  intersect  the  traces  at  a  constant  height  (representing  the 
target  height,  or  trace  endpoint  height).  The  p(ao)  values  at  the  intersections  were 
plotted  against  the  oco  values  as  shown  in  the  accompanying  figures. 

In  Figure  4.38  and  Figure  4.39  where  there  is  no  ray  crossover,  the  P  vs.  Oo  plots 
are  continuous.  Note  in  Figure  4.39  the  rightmost  point  is  estimated  since  there  is  no 
trace  in  (a)  that  just  touches  the  hypothetical  target  height.  In  Figure  4.39  there  are  two  P 
values  for  every  clq,  but  the  curve  is  still  continuous.  In  both  of  these  cases,  8p/5oco  is 
defined  for  all  P  so  that  the  estimation  algorithm  can  eventually  iterate  to  the  correct 
takeoff  angle.  For  instance,  using  Figure  4.39,  suppose  the  target  subtense,  px,  is  0.87 
degrees.  Suppose  further  that  the  initial  takeoff  angle  estimate,  oco,  is  -0.325  degrees 
corresponding  to  two  target  height  (estimated)  subtenses,  pe,  of  0.7  and  0. 18.  Using  the 
solution  that  puts  pe  nearer  Px,  the  slope  of  the  curve  tells  us  to  decrease  oco  to  get  closer 
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to  the  target  subtense.  Continuing  in  this  manner,  specifically  using  equations  3-45,  3-50, 
and  3-51,  the  correct  takeoff  angle  eventually  will  be  found. 


(b) 


Figure  4.38  Subtense  vs.  Takeoff  Angle  Analysis  (Radar  Low,  Target  High,  No 
Ducting),  (a)  traces  for  various  takeoff  angles,  (b)  ray  subtense  at  4000  feet  for  various 
takeoff  angles 
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Subtense  at  Target  Height,  beta  (degrees) 


Figure  4.39  Subtense  vs.  Takeoff  Angle  Analysis  (Radar  High,  Target  Low,  No 
Ducting),  (a)  traces  for  various  takeoff  angles,  (b)  ray  subtense  at  700  feet  for  various 
takeoff  angles 
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Figure  4.40  Subtense  vs.  Takeoff  Angle  Analysis  (Radar  High,  Target  Higher,  Elevated 
Ducting),  (a)  traces  for  various  takeoff  angles,  (b)  ray  subtense  at  3050  feet  (circle)  and 
4500  feet  (star)  for  various  takeoff  angles 


4-38 


Figure  4.41  Subtense  vs.  Takeoff  Angle  Analysis  (Radar  High,  Target  Higher,  Elevated 
Ducting),  (a)  traces  for  various  takeoff  angles,  (b)  ray  subtense  at  350  feet  for  various 
takeoff  angles 
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In  contrast,  the  situations  depicted  in  Figure  4.40  and  Figure  4.41  cause  problems 
for  the  estimator.  Figure  4.40(b)  is  an  example  of  |3  vs.  cxo  curves  with  discontinuities 
caused  by  the  radio  hole.  By  themselves,  discontinuities  are  not  a  problem.  They  make 
certain  subtenses  impossible  to  reach,  but  that  is  a  natural  situation  we  want  to  model. 
However,  the  slope  reversal  at  oto  « -0.35,  caused  by  the  crossover  effects  of  the  duct 
causes  the  estimator  to  go  in  the  wrong  direction.  For  example,  suppose  the  correct 
takeoff  angle,  oco ,  is  -0.1  corresponding  to  a  target  subtense,  Pt  ,  of  1.13  degrees  (in  the 
4500  feet  case),  and  that  the  initial  estimate  is  -0.4  degrees  resulting  in  an  initial 
estimated  subtense,  pe,  of  2.08  degrees.  In  trying  to  following  the  slope  of  the  curve 
down  to  the  target  subtense,  the  estimator  will  eventually  bottom  out  at  some  minimum 
subtense  around  2. 1  degrees.  Of  course,  if  the  initial  guess  is  chosen  close  enough  to  the 
final  takeoff  angle,  the  estimation  will  converge  nicely.  As  written,  however, 
CLIMAREF  estimates  based  on  four-thirds  earth  exclusively. 

In  the  case  of  ducting  the  P  vs.  oto  plot  is  even  more  contorted.  Each  pair  of  lines 
(one  line  in  negative  halfplane  and  one  line  in  positive  halfplane)  represents  one  set  of 
intersections  of  the  vertically  oscillating  rays  with  the  hypothetical  target  height. 
Intuitively,  given  a  small  enough  takeoff  angle,  the  curves  for  each  pair  of  adjacent  cyeles 
will  join  as  approximated  by  the  dashed  lines.  The  single  common  point  of  intersection 
occurs  when  the  trace  just  barely  touches  the  target  height  (instead  of  crossing  it  twice)  as 
illustrated  by  the  dashed  hypothetical  trace  in  Figure  4.41a. 

Examination  of  Figure  4.41b  reveals  four  peculiar  characteristics  of  the  ducting 
situation.  First,  traces  with  small  enough  takeoff  angles  will  never  reach  the  target  height 
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at  all.  Second,  a  given  subtense  may  be  reached  by  traces  of  more  than  one  takeoff  angle. 
Third,  a  given  takeoff  angle  can  reach  multiple  subtenses  (similar  to  the  long,  downward 
path  illustrated  in  Figure  4.39).  And,  finally,  the  curves  are  discontinuous,  not  in  such  a 
way  as  to  absolutely  forbid  reachability  of  certain  subtenses,  but  to  make  the  reachability 
of  some  subtenses  dependent  upon  the  initial  oco  estimate.  For  example,  if  an  initial  Oo 
estimate  of  0.1  is  chosen  and  the  target  subtense  is  0.75,  there  is  no  curve  (no  5p/5ao)  the 
estimator  can  use  to  arrive  at  the  necessary  Oo-  This  is  because  arbitrarily  large  positive 
or  negative  angles  will  not  oscillate.  That  is,  large  positive  angles  leave  the  duct  and 
large  negative  angles  hit  the  ground  to  be  absorbed  or  reflected.  Therefore  the  curves 
ultimately  break  off  on  either  side  making  it  impossible  for  some  initial  takeoff  angles  to 
intersect  a  cycle-pair  that  will  cover  a  subtense  range  containing  the  target  subtense. 
CLIMAREF  is  designed  to  accommodate  the  first  three  difficulties.  This  last  one, 
however,  it  can  not  overcome  using  only  the  four-thirds  earth  estimation. 

So  CLIMAREF  was  tested  for  algorithm  validity,  improvement  over  standard 
atmosphere  and  the  ability  to  trace  in  a  variety  of  scenarios.  Chapter  5  contains 
conclusions  drawn  from  these  results. 
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V.  Findings  and  Conclusions 


Attempting  to  simulate  nature  on  a  digital  computer  is  a  risky  project  comparable 
in  its  audacity  to  building  the  Tower  of  Babel.  The  test  results  described  in  the  previous 
chapter  demonstrate  that  CLIMAREF,  the  model  developed  for  this  thesis,  can 
successfully  predict  radar  beam  paths  and  height  error  with  more  accuracy  than  four- 
thirds  earth  and  standard  atmosphere.  Additionally,  they  indicate  areas  of  the  modeling 
problem  that  require  more  study.  At  this  point,  an  attempt  is  made  to  draw  some 
conclusions  from  those  results  in  order  to  guide  those  whose  task  it  may  be  to  complete 
the  Tower. 

Validation:  CLIMAREF  vs.  Blake  and  CLIMAREF  vs.  RAYS 

The  test  against  Blake  is  particularly  significant  since  NAIC  from  the  beginning 
designated  Blake  as  the  standard  reference.  Happily,  the  simple  test  passed  with  flying 
colors  proving  that  the  algorithms  in  CLIMAREF,  at  least  for  simple  situations,  are  as 
good  as  the  best. 

Likewise,  the  close  match  between  RAYSD  and  RAYS  shows  that  CLIMAREF 
matches  the  industry  standard  as  well  as  can  be  determined  given  the  limited  output 
capability  of  RAYS. 
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Minimum  Elevation  Angle 

The  tests  deseribed  in  Chapter  4  indicate  CLMAREF  will  not  be  able  to 
complete  a  trace  if  the  initial  takeoff  angle  is  so  small  that  the  N-gradient  causes  it  to 
reach  an  extrema  in  the  first  layer  in  which  bending  is  calculated.  Although  time  did  not 
permit  extensive  testing  to  determine  the  relationship  between  takeoff  angles,  layer 
thicknesses,  and  N-gradients,  a  few  general  conclusions  can  be  made. 

The  takeoff  angle  used  in  the  small  negative  angle  test  was  ao=  -5.6530x10'^. 
CLIMAREF  chose  this  angle  to  launch  a  downward  trace  that  would  pass  through  a 
tangent  point  and  hit  the  target  on  the  upswing.  The  target  was  five  meters  above  the 
radar  at  a  range  of  ten  kilometers.  When  the  target  was  raised  one  meter  or  lowered  five 
meters,  the  problem  went  away.  Although  the  specified  target  height  and  range  do  not 
directly  correspond  to  the  initial  estimated  takeoff  angle  (since  this  angle  only  defines  the 
first  of  several  iterations  needed  to  find  the  correct  angle),  some  conclusions  can  be 
drawn.  For  instance,  it  can  safely  be  said  that  the  small  angle  problem  at  a  range  of  ten 
kilometers  is  only  a  problem  when  the  target  height  is  within  a  few  meters  of  the  radar 
height.  Extrapolating  to  a  range  of  one  hundred  kilometers,  we  might  expect  small  angle 
problems  within  a  few  tens  of  meters  of  the  radar  height.  Obviously,  the  significance  of 
the  problem  depends  on  the  application. 

Although  this  limited  test  found  no  problems  resulting  from  a  small  positive 
angle,  the  same  effects  described  in  Chapter  4  for  small  negative  angles  will  presumably 
be  experienced  using  small  positive  angles  within  a  duct.  This  is  because  the  steep 
negative  M-gradients  in  the  upper  portion  of  a  duct  will  cause  a  positive-going  ray  to 
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reach  an  extrema  and  bend  downwards  -  behaving  similarly  to  a  downward  path  in  a  non- 
ducted  atmosphere. 

Additional  measures  need  to  be  devised  to  allow  CLMAREF  to  handle  these 
small  angles.  A  simple,  if  costly,  solution  is  to  decrease  the  layer  thickness  to  something 
less  than  one  meter  (as  it  is  currently  set).  This  reduction  of  the  layer  thickness  will  allow 
smaller  angles  to  be  used,  but  it  will  not  eliminate  the  problem  altogether.  Furthermore, 
the  computing  time  will  increase  by  approximately  the  layer  reduction  factor  since  that 
many  more  layers  will  have  to  be  traced  through  to  reach  a  solution. 

A  similar,  but  more  efficient  solution  might  be  to  reduce  the  layer  thickness  on 
the  fly  only  when  necessary  to  avoid  the  small  angle  problem.  It  may  take  more  iteration 
to  find  the  necessary  layer  thinness,  but  extra  time  only  would  be  required  when  a 
problem  exists.  For  the  vast  majority  of  raytracing  situations  there  will  be  no  small  angle 
problem. 

A  third  solution,  simpler  but  less  effective,  is  to  artificially  eliminate  the  problem 
by  disallowing  small  angles.  Anytime  a  small  takeoff  angle  coming  out  of  the  estimator 
causes  a  problem,  it  could  be  amplified  just  enough  to  eliminate  the  error.  The  negative 
impact  is  that  such  a  solution  would  likely  create  an  artificial  radio  hole,  making  some 
targets  impossible  to  see. 

In  any  case,  small  angles  are  currently  a  weak  spot  in  CLIMAREF  and  more 
testing  will  need  to  be  performed  to  determine  how  best  to  solve  the  problem. 
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Climate  Variation 


The  climate  variation  test  proves  without  a  doubt  that  climatology-based 
raytracing  outperforms  the  standard  atmosphere  model  with  respect  to  accuracy  of 
prediction.  With  the  target  at  10,000  feet  and  60  nautical  miles  out,  CLIMAREF 
predicted  height  error  magnitudes  ranging  from  101.2  to  2256.0  feet.  These  are, 
respectively,  702.8  and  1452.0  feet  away  from  the  height  error  calculated  using  standard 
atmosphere,  804.0  feet.  Although  the  overall  mean  of  987.0  feet  was  close  to  standard 
atmosphere,  overall  average  deviation  from  standard  atmosphere  was  about  400  feet.  A 
particularly  unusual  result  was  calculated  for  Fort  Lamy,  Chad  with  elevated  ducting  in 
January:  the  height  error  was  actually  negative  at  -1973. 1  feet.  Such  radical  departures 
from  standard  atmosphere  will  be  useful  to  radar  engineers  using  AMBER  to  more 
accurately  predict  radar  performance. 

The  climate  variation  test  provides  insight  not  only  into  the  utility  of 
CLIMAREF,  but  also  into  the  propagation  effects  that  can  be  expected  under  various 
climatic  conditions.  Note  that  these  results,  as  children  of  the  HEPC  database,  represent 
mean  conditions  only.  No  variation  data  is  available  or  represented. 

In  addition  to  proving  the  superior  accuracy  of  CLIMAREF,  the  data  provides 
many  insights  into  the  climatological  aspects  of  this  subject.  For  example,  the  extreme 
climates  of  Jodhpur  (monsoon  region),  Chad  (the  Sudan),  and  the  coast  of  Libya 
(Mediterranean  climate)  are  reflected  in  the  wide-ranging  height  errors.  Particularly 
notable  is  the  negative  height  error  corresponding  to  January  in  Chad  under  elevated 
ducting  conditions.  This  unusual  effect  occurred  because  the  ray  was  subrefracted,  that 
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is,  it  actually  bent  away  from  the  earth.  Obviously,  the  superrefractive  case  is  more 
typical,  but  subrefraction  can  and  does  occasionally  occur.  On  the  other  hand,  places  with 
more  moderate  climates  like  Port  Hardy,  Churchill,  Denver,  Bismarck,  and  Columbia  are 
notable  for  their  lack  of  variation  and  for  how  closely  they  agree  with  the  standard 
atmosphere. 

The  conclusion  is  obvious.  Climate  does  affect  radar  beam  bending  to  a  large 
extent.  An  appreciation  for  how  significant  the  effect  is  may  be  gained  by  considering  the 
Federal  Aviation  Administration  mandates  1000  feet  vertical  separation  between  air 
traffic  patterns  as  an  acceptable  safety  margin.  The  variation  of  the  error  data  presented, 
being  on  the  same  order  of  magnitude,  makes  the  effects  of  climate  worth  noting.  Of 
course,  the  significance  ultimately  depends  on  the  application  of  the  model,  but  the  FAA 
standard  is  a  good  benchmark. 

Capability  of  Model  to  Handle  All  Types  of  Refractive  Effects 

These  tests  prove  that  CLIMAREF  can  model  most  radar-target-atmosphere 
scenarios  with  only  a  few  exceptions.  Again,  the  model  can  recognize  when  the  target  is 
in  a  radio  hole  and  when  it  is  over  the  horizon.  It  can  trace  upwards  or  downwards  and 
can  properly  work  through  the  tricky  business  of  a  tangent  point.  Furthermore,  it  can 
differentiate  between  a  target  on  the  short  downward  path  or  the  long  downward  path. 
Also,  it  can  model  many  ducting  effects,  albeit  with  a  few  more  iterations  of  the 
estimator.  It  works  properly  for  the  most  common  situations. 
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Nonetheless,  CLIMAREF  is  not  perfect.  It  can  not  satisfactorily  trace  when  the 
duct  causes  the  traces  to  cross  or  when  the  path  to  the  target  is  a  pure  ducted  path. 

Probably  the  most  striking  observation  about  both  situations  is  that  the  second 
limitation  of  raytracing  (discussed  in  Chapter  2)  is  being  violated:  that  is,  neighboring 
rays  must  remain  close  to  parallel  within  one  wavelength.  Further  research  is  required  to 
determine  exactly  why  the  model  atmosphere  violates  this  condition.  Most  likely,  the 
steep  gradients  inherent  in  the  inversions  are  the  problem.  Depending  on  the  frequency  in 
use,  these  gradients  may  also  violate  the  first  condition  of  geometric  optics  (the  refractive 
index  must  not  vary  appreciably  in  one  wavelength).  So,  it  may  be  that  completely 
different  methods  need  to  be  used  to  model  the  more  serious  effects  of  ducting. 

On  the  other  side,  EREPS  RAYS,  which  is  a  component  of  the  premiere 
modeling  package  for  this  sort  of  thing,  traces  in  and  through  the  ducts  matching 
CLIMAREF  trace  for  trace,  with  rays  oscillating  and  crisscrossing  all  through  the  duct. 

So,  either  both  models  are  wrong  to  raytrace  in  the  duct,  or  the  limitations  do  not  apply  in 
this  case  (for  some  reason)  and  problems  in  CLIMAREF  must  be  overcome  another  way. 

Apart  from  considerations  of  the  applicability  of  raytracing,  the  P  vs.  oco  plots  at 
the  end  of  the  last  chapter  suggest  that  part  of  the  problem  has  to  do  with  the  initial 
elevation  angle  estimate.  Both  scenarios  in  which  tracing  will  potentially  break  down 
(Figure  4.40,  and  Figure  4.41)  can  be  made  to  work  if  the  initial  takeoff  angle  estimate  is 
close  enough  to  the  final  required  value.  By  taking  the  effects  of  the  ducts  into  account,  it 
is  possible  a  more  sophisticated  estimation  procedure  can  be  developed  that  will  make 
these  disturbing  anomalies  irrelevant.  Alternatively,  logic  could  be  employed  to 
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recognize  the  traps  caused  by  such  effects  in  order  to  redirect  the  estimation  and  iteration 
procedure  around  them. 

Nevertheless,  the  question  of  whether  or  not  these  techniques  violate  the  basic 
assumptions  of  geometric  optics  must  be  answered  and  the  answers  must  be  dealt  with. 
Since  most  analysis  of  ducting  is  performed  using  waveguide  theory  and  physical  optics, 
it  may  be  that  these  disciplines  are  the  only  way  to  find  the  path  of  the  energy  through  the 
duct,  if  indeed  a  well-defined  path  exists  at  all.  While  CLIMAREF  as  it  stands  is  a 
worthy  and  useful  tool  because  of  the  wide-variety  of  scenarios  it  can  model,  further 
study  will  be  required  to  resolve  these  questions  if  it  is  ever  to  be  useful  and  reliable 
under  all  circumstances. 

Recommendations  for  Further  Research 

Due  to  time  limitations,  CLIMAREF  was  not  developed  as  fully  as  it  could  be. 
No  doubt  many  improvements  can  be  made.  Here  are  a  few  ideas. 

First,  compare  CLIMAREF’ s  predictions  against  actual  height  error  statistics. 

The  great  difficulty  in  doing  this  ultimate  test  is  getting  the  data.  Most  ground-based 
radars  already  make  an  attempt  to  compensate  for  beam  bending.  Therefore  the  height- 
measurements  collected  and  stored  are  not  pure.  Possibly,  the  raw  data  could  be  tapped 
from  the  radar  receiver  during  a  dedicated  test.  If  enough  data  were  taken,  statistics  could 
be  derived  and  compared  to  CLIMAREF  predictions.  Such  a  test  would  be  a  true 
measure  of  the  model’s  worth. 
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Deciding  how  to  trace  through  both  the  surface  duct  and  the  elevated  duct  at  the 
same  time  could  enhance  the  worth  of  the  model.  First,  the  double  duct  will  have  to  be 
properly  constructed.  Patterson  (1987:  15-16)  describes  the  construction  of  each  duct 
separately,  but  says  nothing  about  two.  Then,  all  the  issues  regarding  radio  holes  and  the 
unique  paths  produced  by  such  a  combination  will  have  to  be  studied  and  accommodated 
with  additional  logic  built  into  the  model. 

Another  idea  is  to  apply  the  interpolation  routines  to  the  ducting  statistics  as  well 
as  the  surface  statistics.  Care  must  be  taken  to  determine  when  these  statistics  can  and 
can  not  be  averaged. 

To  make  the  climatology  database  more  accurate  and  useful,  several  tasks  might 
be  accomplished:  Variation  statistics  could  be  researched  and  added  to  the  mean  data  in 
the  database.  This  would  give  the  modeler  a  better  idea  of  the  range  of  height  error  he 
can  expect  in  a  given  climatic  situation.  Also,  the  database  could  be  expanded  to  provide 
a  finer  resolution  of  radiosonde  stations  in  certain  areas  of  interest. 

Finally,  to  be  really  useful,  CLIMAREF  could  take  into  account  reflections  and 
multipath  and  compute  the  power  loss  to  the  target  using  physical  optics. 

Conclusion 

CLIMAREF,  then,  has  been  found  to  be  a  good  tool  for  simulating  radar  beam 
bending.  It  performs  well  in  most  situations,  and  breaks  down  in  some.  However,  the 
troubling  issues  do  not  seem  to  be  irresolvable.  This  exploration  of  the  issues  involved 
with  climatology-based,  geometric  optics-based  prediction  of  radar  beam  bending  is  by 
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no  means  the  first  or  the  last  word  on  the  subject.  However,  it  is  the  author’s  hope  that  it 
will  play  a  small  part  in  advancing  the  state-of-the-art  of  radar  modeling. 
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Appendix  A 


CLIMAREF  MATLAB  Code 

Note:  Because  of  wraparound  some  of  this  code,  as  printed,  is  improperly  formatted 

function  CRlmain 

global  MON  DYNT  RLAT  RLONG  REL  RAZ  THT  TRG 

global  ISOLATED  USESTAT 

global  NS  NIK  SEL  DUCTFLG  MULTSTATFLG 

global  GMSB  OHSB  MTSB  MDSB  MFSB  PSB  GMEL  OREL  MTEL 

global  MDEL  MFEL  PEL  P2EL  PSBEL  PHs  PMs  PHe  PMe  PHn  PMn 

global  RO  HEIGHTS  BETAES  TRGAPP  ALFG  ALFO 

global  THTAPP  THTERR  THTERRIOO  TRGERR  TRGERRIOO  DUCTHTS  ERRMSG 
%CLIMAREF  Main  module,  "crlmain.m" 

%This  is  the  main  loop.  It  controls  all  other  processes 
% 

%VARIABLES : 

%ALF0  -  takeoff  angle  of  refracted  path  to  radar 
%ALFG  *-  geometric  elevation  angle  of  target  at  radar 
%ang  -  multipurpose  angle  variable 

%another  -  user  input  whether  or  not  to  run  program  again  (string:  y/n) 
%beta  -  subtense  counter 

%BETAES  -  list  of  betas  defining  subtenses  of  every  point  on  the  trace 

%betastep  -  rgstep  in  terms  of  subtense 

%DUCTFLG  -  indicates  type  of  ducting  user  chose 

%endloop  -  flag  to  signify  user  termination  of  the  program 

%ERRMSG  -  error  message  (if  any)  output  from  CR9raytrace 

%HEIGHTS  -  list  of  heights  of  every  point  on  the  raytrace 

%htsft  -  htskm  converted  to  feet 

%htskm  -  heights  in  kilometers,  used  for  plotting  geometric  path  on  flat 
earth  plot 

%htstep  -  distance  between  height  grid  lines  on  curved  earth  plot 
%intype  -  input  data  by  file  or  manually 

%ISOLATED  -  flag  indicates  radar  site  has  no  radiosonde  stations  around 
it 

%maxhtft  -  maximum  height  from  surface  ray  reaches  (in  feet) 

%maxrad  -  maximum  distance  from  center  of  earth  that  the  ray  reaches 

%maxrng  -  maximum  ground  range  of  target 

%maxrngNmi  -  maximum  ground  range  of  target  in  Nmi 

%MULTSTATFLG  -  indicates  whether  user  chose  2  or  3  stations  to 

interpolate 

%N1K  -  refractivity  at  1km  above  surface 

%NS  -  refractivity  at  surface 

%p,q  -  generic  cartesian  coordinates 

%PEL  -  probability  of  an  elevated  duct  present  (from  duct  stats) 
%plottype  -  flat  earth,  *f',  or  curved  earth,  ’c* 

%PSB  -  probability  of  a  surface  duct  present  (from  duct  stats) 

%R0  -  radius  from  center  of  earth  to  surface  (=a+SEL) 

%REL  radar  elevation  in  kilometers 
%relft  -  radar  elevation  in  feet 

%rgstep  -  distance  (on  curved  earth  plot)  between  ground  range  ticks 
%seeplot  -  y  or  n,  whether  or  not  user  wants  to  see  the  plot 
%tabeta  -  subtense  to  target  based  on  apparent  range  and  ray 
%  initial  angle  of  elevation 

%tbeta  -  subtense  to  target  based  on  geometric  range  and  geometric 
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%  initial  angle  of  elevation 

%THTAPP  -  apparent  (measured)  target  height 
%ticklength  -  length  of  ground  range  tick  marks 
%TRG  -  actual  target  range 

%TRGAPP  -  apparent  (measured)  target  range 
%trgparts  -  geometric  range  segments 

%useduct  -  s,e,  or  n^  whether  or  not  user  wants  to  use  duct  info  for 
%  the  raytrace 

%x,  y, xl, yl ^ x2, y2  -  multipurpose  cartesian  coordinates 

Q. 

■& 

% 


% 

Q, 

'O 

endloop=0; 
another= [ ] ; 
intype=[] ; 

MULTSTATFLG=0; 
disp ( '  * ) 

while  isempty (intype) 

intype=input ( *  Input  parameters  Manually  or  from  the  preset  File 
(m/f) :  ’,^s'); 

if  isempty  ( intype )  |  (intype'-=*m’ &intype'-=*M^  &intype'-=*  f  ^ &intype~=' FM 
disp(’M  or  F,  please.*) 
intype=  [] ; 

end 

end 

if  intype==' F*  I intype==* f * 

CR2in_pre  %Call  pre-set  user  input  module 
else 

CR2in  %Call  user  input  module 

end 

CR3f indnerstat s  %Find  nearest  stations  to  radar 

if  ISOLATED==0  %Run  these  if  the  radar  DOES  have  stations  around  it 
CR4calcdist  %Find  distances  of  nearest  stations 
CRSpickstat  %Allow  user  to  choose  a  radiosonde  station 
CRGloadNdata  %Retrieve  ref ractivity,  elevation  values  for  station 

else 

disp ( *  * ) 

disp (* Your  radar  is  too  far  away  from  any  radiosonde  station  in  the') 
disp (' database .  The  refractivity  profile  will  be  calculated  using') 
disp ('the  standard  atmosphere  (universal  average).') 

NS=313; 

NIK-271; 

end 

if  length (NS) >1  %If  loadNdata  returned  more  than  one  NS, 

MULTSTATFLG-1; 

CRVinterp  %interpolate 

end 

if  NS—0 

disp ( '  ' ) 

disp ('No  data  for  this  station.  Run  again  and  choose  a  different 
station. ' ) 

else 
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CRSductstats  %get  duct  statistics 

if  -PSB&~PEL, disp ( ’No  ducting  data  for  this  station end 
if  USESTAT'-=999&-'MULTSTATFLG&  (PSBI  PEL) 
disp ( ’  ’ ) 

disp( ’SURFACE-BASED  DUCTS’) 

disp  ([’ Surf  ace  to  Inflection  Pt  N-unit  Gradient:  ’  nuin2str  (GMSB)  ]  ) 

disp ([ ’Optimum  Ht  (layer  base):  ’  num2str (OHSB)  ’km’]) 

disp ([’Duct  Thickness:  ’  num2str (MTSB)  ’  km’]) 

disp ( [ *M-Unit  Deficit:  '  num2str (MDSB) ] ) 

disp ([’Max.  Frequency  Trapped:  ’  num2str (MFSB)  ’  MHz’]) 

disp ([* Percent  of  Time  Duct  Occurs:  ’  num2str(PSB)  ’  %’]) 

disp ( ’  ’ ) 

disp ( ’ELEVATED  DUCTS’) 

disp ([* Surface  to  Inflection  Pt  N-unit  Gradient:  ’  num2str (GMEL) ] ) 
disp ([ ’Optimum  Ht  (layer  base):  '  num2str (OHEL)  ’  km’]) 
disp([’Duct  Thickness:  ’  num2str (MTEL)  ’  km’]) 
disp ( [ ’ N-Unit  Deficit:  ’  num2str (MDEL) ] ) 

disp([’Max.  Frequency  Trapped  (MHz):  ’  num2str (MFEL)  ’  MHz’]) 
disp ([’ Percent  of  Time  Duct  Occurs:  ’  num2str(PEL)  ’%’]) 
disp ([’ Probability  of  >1  Elevated  Duct:  ’  num2str (P2EL)  ’  %’]) 

^  disp ( ’  ' ) 

disp ([’ Probability  of  Surface-Based  and  Elevated  Duct:  ’ 
num2str ( PSBEL)  ’  %’]) 

useduct= [ ] ; 
disp ( ’  ’ ) 

while  isempty (useduct ) 

disp (’Take  into  account  ducting:’) 
if  PSB&PEL 

useduct=input (’ Surface  duct,  Elevated  duct,  or  No  duct 
(s/e/n) ?  ’ , ’s’ ) ; 

if 

isempty  (useduct )  |  (useduct~=’ s  ’  Suseduct'-^’  S  '  &useduct~=’  e*  &  .  .  . 

useduct~= ’ E ’ &useduct^- ’ n ’ &useduct~= ’ N ’ ) 
disp(’S,  E,  or  N,  please.’) 
useduct= [ ] ; 

end 

elseif  -PSB&PEL 

useduct=input (’ Elevated  duct  or  No  duct  (e/n)?  ',’s’); 
if  isempty  (useduct )  |  (useduct ’  e  ’  &useduct~=  ’  E .  .  . 
useduct~=  ’  n '  &useduct'-=  ’  N  ’  ) 
disp('E  or  N,  please.*) 
useduct= [ ] ; 

end 

elseif  PSB&'-PEL 

useduct=input (’ Surface  duct  or  No  duct  (s/n)?  ’,'s'); 
if  isempty (useduct) | (useduct~= ’ s ’ &useduct-= ’ S ’ & . . . 
useduct-'=  ’  n  ’  &useduct'-=  ’  N ' ) 
disp('S  or  N,  please.’) 
useduct^ [ ] ; 

end 

end 

end 

if  useduct==’S’ | useduct== ’ s ’ 

DUCTFLG=1; 
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elseif  useduct==’ E ' 1 useduct==’ e ' 

DUCTFLG=2; 

elseif  useduct== ’ N ’ I useduct== * n * 

DUCTFLG=0; 

end 

else 

DUCTFLG^O; 

end 

CR9raytrace  %do  raytracing  to  construct  ray  path 

%Display  Calculated  Values  to  MatLab  Command  Window 
disp ( ’  ’ ) 

disp(*  ’) 

disp ([ ’Actual  Height:  ’  num2str (THT*3280 . 84  )  ’ft’]) 

disp {[ ’Apparent  Height:  ’  num2str (THTAPP*3280 . 84 )  'ft’]) 

disp( [’Height  Error:  ’  num2str (THTERR*3280 . 84 )  'ft']) 

disp ([’ Height  Percent  Error:  ’  num2str (THTERRIOO ) ] ) 
disp(’  ’) 

disp  ([ ’Actual  Range:  '  num2str  (TRG"^  .  53996)  'Nmi’]) 

disp ([ ’Apparent  Range:  ’  num2str (TRGAPF* . 53996)  ’Nmi’]) 

disp ([’Range  Error:  ’  num2str (TRGERR* . 53996)  ’Nmi’]) 

disp( [’Range  Percent  Error:  ’  num2str (TRGERRIOO) ] ) 

%Plot  and  Display  Information 

if  isempty (ERRMSG) 
plottype=[] ; 
disp { '  ’ ) 

while  isempty (plottype) 

plottype=input (’ Curved  earth,  Flat  earth,  No  plot  (c/f/n) ? 

^  *  s  ’  )  ; 

if 

isempty  (plottype)  |  (plottype-=  ’  c  ’  Scplottype-^  ’  C  ’  &plottype~=  ’  f  ’  &plottype~=  ’ 
F  ’  &plottype'-= '  n  '  &plottype-'= '  N  ’  ) 

disp(’C,  F,  or  N,  please.’) 
plottype= [] ; 

end 

end 

else  %if  there’s  an  error  message,  don’t  plot  anything 
plottype= ’ n ’ ; 
disp ( ’  ’ ) 

disp (’!!!!!!!!!!  i  1  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'  ) 
disp (ERRMSG) 

disp ('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  1  !!!!!!!!!!!!'  ) 
disp  (  '  ’ ) 

end 


relft=REL*3280. 84; 

if  plottype===  ’  C  ’  I  plottype==  ’  c  ’ 

%plot  bent  path 
ang-[ (pi/2 ) -BETAES ] ; 
x= [RO+HEIGHTS] . *cos (ang) ; 
y-[RO+HEIGHTS]  .  ’^sin(ang)  ; 
plot (x,y, ’r’ ) 
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hold  on 


%plot  geometric  path 
p=[x{l) ,x(l)+TRG*cos (ALFG) ] ; 
q=[y(l) ,y(l)+TRG*sin(ALFG) ] ; 
plot (p,q, ’c’ ) 

%plot  straightline  path  based  on  TRGAPP  and  inital  ray  angle, ALFO 
p=[x(l) ,x{l)+TRGAPP*cos (ALFO)  ]  ; 
q=[y(l) ,y{l)+TRGAPP* sin (ALFO) ] ; 
plot (p,q, ) 

%plot  earth  surface  (technically,  the  local  geoid  —  using  radius 

%of  curvature  instead  of  the  actual  local  earth’s  radius) 

ang=[pi/2-BETAES (length (BETAES) ) : .0001:pi/2] ; 

x=R0*cos (ang) ; 

y==R0^sin  (ang)  ; 

plot (x, y) 

%plot  concentric  circles  around  earth 
ang-[pi/2-BETAES ( length (BETAES ) ) : .0001:pi/2] ; 
if  THTAPP>REL 

maxrad=RO+THTAPP; 

else 

maxrad“RO+REL; 

end 

maxht ft= (maxrad-RO ) *3280.84; 
if  maxhtft>=30000 
htstep==10000; 

el seif  maxhtft<30000&maxhtft>=10000 
htstep=5000; 

el seif  maxhtft<10000&raaxhtft>=5000 
htstep=1000; 
else 

htstep=500; 

end 

for  r= [RO : htstep/3280 . 84 imaxrad] 
x=r*cos (ang) ; 
y=r*sin (ang) ; 
plot (x, y, ’b: ’ ) 

end 


%Plot  ground  range  tick  marks  &  labels  on  earth  surface 
maxrng=TRG*cos (ALFO) ; 
maxrngNmi=maxrng* .53996; 

if  maxrngNmi>30 
rgstep=10 ; 

el seif  maxrngNmi<=30&maxrngNmi>15 
rgstep=5 ; 
else 

rgstep=l; 

end 

betastep=rgstep/ (RO* . 53996)  ; 
ticklength^. 02* (maxrad-RO)  ; 
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for  beta= [0 :betastep:maxrng/RO] 
xl==R0*cos  (pi/2-beta)  ; 
yl=R0*sin  (pi/2“-beta)  ; 
x2= (RO+ticklength) *cos (pi/2-beta) ; 
y2= (RO+ticklength) *sin (pi/ 2 -bet a) ; 
plot ( [xl,x2] , [yl,y2] ) 
if  beta 

text (x2, y2, num2str (beta^RO^ . 53996) , ' VerticalAlignment ' , . . . 
’bottom' HorizontalAlignment Center M  %add  Nmi  labels 

end 

end 

%plot  duct  boundaries  if  applicable 
if  DUCTFLG 

ang=[pi/2-BETAES (length (BETAES) ) : .0001:pi/2] ; 
for  ductht= [ 1 : 2 : 3] 

if  ductht  %i.e.  don't  plot  a  zero-level  duct  boundary 
x=(DUCTHTS (ductht) +R0) ^cos (ang) ; 
y- (DUCTHTS (ductht ) +R0 ) ^sin (ang) ; 
plot  (X,  y,  'g— '  ) 

end 

end 

end 

%arrange  and  label  axes 
axis  tight 

set (gca,  ' ytick  * ,  [RO : htstep/3280 . 8  4 : maxrad] ) 

set (gca, ' ytickLabel ' , [0 :htstep:maxhtft] ) 

set (gca, ' xtick' , [ ] ) 

xlabel (' Ground  Range  (Nmi)') 

ylabel ( 'Altitude  (feet)') 

set (gca, ' FontSize ' , 8 ) 

set (gca^  ' FontName Times  New  Roman') 
hold  off 

elseif  plottype== ' F' I plottype== ' f '  %Plot  type  is  flat  earth 
%plot  bent  path 

gndrngs=BETAES*RO* . 53996/  %conversion  factor  changes  km  to  Nmi 
htsft=HEIGHTS*3280 . 84 ; 
plot (gndrngs, htsft , 'r' ) 
hold  on 

%plot  geometric  path 
trgparts=[0:TRG/length (BETAES) :TRG] ; 

htskm=  (sqrt  (  (RO+REL)  ^2+trgparts  .  '^2-2’^  (RO+REL)  ^trgparts*  .  .  . 
cos (pi/2+ALFG) )-R0) ; 

gndrngs-RO* . 53996*asin ( trgparts^sin (pi/2+ALFG) . / (RO+htskm) ) ; 
htsf t=htskm*3280 . 84 ; 
plot (gndrngs, htsft , ' - . c ' ) 

%plot  straightline  path  based  on  TRGAPP  and  inital  ray  angle, ALFO 
trgparts=[0:TRGAPP/ (length (BETAES ) -1 ) iTRGAPP] ; 
htskm=  (sqrt  (  (RO+REL)  '^2+trgparts  .  -^2-2*  (RO+REL)  ^trgparts*  .  .  . 
cos (pi/2+ALFO) ) -RO) ; 

gndrngs=R0^ . 53996*asin ( trgparts^sin (pi/2+ALFO ) . / (RO+htskm) ) / 
ht sf t=htskm'*'3280 . 84  ; 
plot (gndrngs , htsft , ' - .m' ) 
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%plot  duct  heights  (if  applicable) 
if  DUCTFLG 

plot ( [ 0 , gndrngs ( length ( gndrngs ))],... 

[DUCTHTS  (1)  *3280. 84,  DUCTHTS  (1)  ^3280. 84]  ,  g’  ) 
plot ( [ 0 , gndrngs ( length ( gndrngs ))],... 

[DUCTHTS (3) *3280.84, DUCTHTS (3 ) *3280 . 84 ] , '--g* ) 

end 

xlabel ( ^ Ground  Range  (Nmi)’) 

ylabel (’ Height  above  sea  level  (ft)’) 

legend ( ’Actual  Path’, 'Direct  Path ', ’Apparent  Path’, 'Duct 
Boundaries ' ) 
end 

end 

clear 
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%CLIMAREF  User  Input  module  ”CR2in.m" 

%This  module  allows  the  user  to  input  the  following  data: 

%month  (integer  1-12) 

%day  or  night  (l=day,  0=night) 

%radar  latitude  in  degrees  (between  -70  and  +80) 

%radar  longitude  in  degrees 

%radar  elevation  entered  in  feet  AGL,  converted  to  km 
%radar  pointing  azimuth  in  degrees  (north=0) 

%actual  target  height  entered  in  feet  AGL,  converted  to  km 

%actual  straight  line  distance,  radar  to  target  in  Nmi  converted  to  km 

Q. 

O 

% INPUT  VARIABLES 
%none 

%  ' 

% INTERNAL  VARIABLES 

%dyntstr  -  string  form  of  DYNT,  can  equal  D,N,  or  B 
%relft  -  radar  elevation  in  feet 
%thtft  -  target  height  in  feet 
%trgnm  -  target  range  in  Nmi 

a 

"D 

%OUTPUT  VARIABLES 

%MON  -  month  of  interest  (integer  1-12) 

%RLAT  -  latitude  of  radar  (-70  to  +80  deg) 

%RLONG  -  longitude  of  radar  (deg) 

%REL  -  elevation  of  radar  (km) 

%RAZ  -  azimuth  radar  is  pointing  to  (deg) 

%THT  -  actual  target  height  (km) 

%TRG  -  actual  target  range  (km) 

redo=l ; 
while  redo==l 


MON=  [  ]  ; 
dyntstr= [ ] ; 

RLAT= [ ] ; 

RLONG- [ ] ; 
relft=[]  ; 

RAZ-  [  ]  ; 
thtft=[]; 
trgnm=[] ; 
crct= [ ] ; 

dispC  ') 

disp(' Please  input  the  following  information:’) 
disp ( '  ’ ) 

while  isempty(MON) 

MON=input (' Month  (1-12):  ’); 

if  isempty(MON) | round (MON) ~=MON | M0N<1 | MON>12 
disp(’You  must  enter  a  month  number  (1-12)’) 

MON=  [  ]  ; 

end 

end 

while  isempty (dyntstr ) 

dyntstr=input ( ’ 1200Z, OOOZ,  or  avg  of  both  (D/N/B) :  ’,’s’); 

if 

isempty  (dyntstr)  |  (dyntstr'-^’  d’  Scdyntstr-^  '  D ’  &dyntstr-=  ’ n  ’  &dyntstr~=  ’ N ’  &dy 
ntstr-'=’b'  &dyntstr~=’ B’  ) 
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disp('You  must  indicate  day, night,  or  both’) 
dyntstr= [ ] ; 

end 

end 

if  dyntstr=='D’ | dyntstr== ’ d ’ 

DYNT=1; 

elseif  dyntstr== ’ N ’ i dyntstr== ' n ’ 

DYNT-0; 

else  %both  day  and  night  (avg) 

DYNT-2 ; 

end 

while  isempty (RLAT) 

RLAT=input {’ Radar  latitude  (-70  to  +80  deg):  ’); 
if  isempty (RLAT) | RLAT<-70 | RLAT>+8 0 

disp('You  must  enter  a  latitude  value  between  -70  and  +80 

degrees ’ ) 

RLAT= [ ] ; 

end 

end 

while  isempty (RLONG) 

RLONG=input (’ Radar  longitude  (+/-  deg):  ’); 
if  isempty (RLONG) | RLONG<-180 | RLONG>180 

disp{’You  must  enter  a  longitude  value  between  -180  and  +180 

degrees ’ ) 

RLONG= [ ] ; 

end 

end 

while  isempty ( relft ) 

relft=input (’ Radar  elevation  (feet  AGL) :  ’); 

if  isempty ( relft ) I relft<0 

disp(’You  must  enter  a  positive  elevation  value’) 
relft= [ ] ; 

end 

end 

REL=relft* . 0003048;  %convert  feet  to  kilometers 
while  isempty (RAZ) 

RAZ=input (’ Radar  azimuth  (deg)  -  default=0  (north):  ’); 
if  isempty (RAZ) 

RAZ=0; 

disp(’The  radar  is  pointing  to  0  degrees  (north)’) 

end 

if  RAZ<0 |RAZ>360 

disp('You  must  enter  an  azimuth  value  between  0  and  360 

degrees ' ) 

RAZ= [ ] ; 

end 

end 

while  isempty (tht ft ) 

thtft=input ( ’Actual  Target  Height  (feet  AGL):  '); 
if  isempty (thtft) I thtft<0 

disp('You  must  enter  a  positive  target  height’) 
thtft= [ ] / 

end 
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end 

THT=thtft* . 0003048 ;  %convert  feet  to  kilometers 
while  isempty ( trgnm) 

trgnm==input  {  ^Actual  Target  Range  (Nmi):  *); 
if  isempty (trgnm) | trgnm<0 

disp(*You  must  enter  a  positive  target  range  value’) 
trgnm=[]  ; 

end 

end 

TRG=trgnm^l . 852 ;  %convert  Nmi  to  kilometers 

%Check  to  see  if  all  the  entries  are  correct 
while  isempty (crct) 
disp(’  ’) 

crct=input ( ’ Are  all  the  above  entries  correct  (y/n) ?  ’,’s'); 
if  isempty  (crct)  |  (crct'-=’ y  ’  &crct~=  ’  Y’ &crct'-=  ’  n '  &crct~=  ’  N '  ) 
disp(’ Please  answer  Y  for  yes,  or  N  for  no.') 
crct= [ ] ; 

end 

end 

if  crct== ' Y ' I crct== ' y ' 
redo=0 ; 

end 


end  %process  will  drop  out  of  while-loop  here  as  long  as  crct  was  'y' 


%Clear  unneeded  internal  variables 
clear  dyntstr  relft  thtft  trgnm 

pack 

%Now  return  to  CRmain 
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%CLIMAREF  Input  module  --  preset  data 

%This  module  allows  the  user  to  input  the  following  data: 

%MON  -  month  (integer  1-12) 

%DYNT  -  day  or  night  (l=day,  0=night) 

%RLAT  -  radar  latitude  in  degrees 
%RLONG  -  radar  longitude  in  degrees 
%REL  -  radar  elevation  in  km  MSL 

%RAZ  -  radar  pointing  azimuth  in  degrees  {north=0) 

%THT  -  actual  target  height  in  km  MSL 

%TRG  -  actual  straight  line  distance  from  radar  to  target  in  Nmi 
% 

M0N=7 ; 

DYNT=0; 

RLAT=42; 

RLONG— 80; 

REL-. 3; 

RAZ=0; 

%THT=3, 048; 

THT=input { *  THT  ( km)  ’ ) ; 

%TRG=111.1; 

TRG=input ( *  TRG  ( km)  ’ ) ; 
disp ( ^  * ) 

disp{^ Radar  Parameters:*) 
disp  (  *  * ) 

disp ([* Month :  '  num2str (MON) ] ) 

disp( [*Day/Night:  *  num2str ( DYNT ) ] ) 

disp (['Radar  Lat :  '  num2str (RLAT) ] ) 

disp (['Radar  Long:  '  num2str (RLONG) ] ) 

disp (['Radar  Elevation  (km):  '  num2str (REL) ] ) 

disp (['Radar  Azimuth:  '  num2str (RAZ ) ] ) 

disp ([' Target  Height  (km):  '  num2str (THT ) ] ) 

disp ([' Target  Range  (km):  '  num2str (TRG) ] ) 
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function  CR3f indnerstats 

global  ISOLATED  RLAT  RLONG  NEARMSQS  NEARSTATS 

%CLIMAREE  Load  climate  data  module,  submodule  1  (CR3f indnerstats . m) 
%This  submodule  find  the  nine  MSQs  surrounding  the  radar  and  the  record 
%numbers  of  all  the  radiosonde  stations  within  those  MSQs.  If  there 
%aren*t  any  it  sets  a  flag. 

% INPUT  VARIABLES: 

%RLAT  -  radar  latitude  (deg) 

%RLONG  -  radar  longitude  (deg) 

%INTERNAL  VARIABLES: 

%lat  -  lat  of  msq  of  interest  (radar  lat  +/-  0  or  10) 

%listind  -  used  to  index  msqlist 

%long  --  long  of  msq  of  interest  (radar  long  +/-  0  or  10) 

%msq  -  msq  number  of  record  currently  being  examined 
%msqind  -  indices  the  nine  msqs  surrounding  the  radar 

%msqlist  -  list  of  all  msqs  in  the  order  they  appear  in  the  data  file, 

%  the  index  of  the  first  record  bearing  the  given  msq  number, 

%  and  the  number  or  records  bearing  that  msq  number  (msqlist 

%  is  stored  as  a  data  file  until  it  is  needed  by  the  program 

%ordflag  -  flag  indicating  near  MSQs  are  in  order 

%pmmsqind  -  index  of  msq  immediately  west  of  PM  from  southernmost 
%  to  northernmost 

%pmmsqs  -  all  the  msqs  immediately  west  of  PM  from  south  to  north 

%sortptr  -  pointer  used  to  track  MSQs  in  bubble  sort 

%swapflag  -  flag  used  in  sort  to  indicate  a  swap  of  positions  has 

happened 

%tempmsq  -  temporary  holder  used  in  sort 
% 

%OUTPUT  VARIABLES: 

%NEARMSQS  -  numbers  of  the  9  marsden  squares  surrounding  the  radar 
%NEARSTATS  -  indices  of  all  stations  in  the  9  marsden  squares 
surrounding 

%  the  radar  (variable  length  vector) 

%ISOLATED  -  flag  indicated  there  are  no  radiosonde  stations  in  any  of 
%  the  nine  MSQs  surrounding  the  radar  (a  value  of  1  indicates 

%  isolation) 

ISOLATED=0; 

9  MSQ's  nearest  the  radar****^***^**-^**^^*^*^^^^-^**^ 
pmmsqs=[516  480  440  408  372  336  300  1  37  73  109  145  181  217  253]; 
msqind=0;  %MSQ  index 

for  lat= [RLAT~10 : 10 : RLAT+10] ;  %cover  all  msq's  surrounding  the 

for  long==  [RLONG-10  : 10  :  RLONG+10  ]  ;  %radar,  3  rows  and  three  columns 
msqind=msqind+l;  %increment  index 

if  long>180  %if  RLONG  near  180  merid,  nearby  msq*s  need  new 
long=long~360 ;  %latitude  reference  calculated 

end 

if  long<-180  %likewise 
long=long+3  60 ; 

end 

if  lat<-70 I lat>80  %if  msq  doesn*t  exist,  set  it  to  0  and 

NEARMSQS (msqind) =0 ;  %go  on  to  the  next  one 
else 

pinmsqind=floor  (lat/10) +8 ;  %find  the  msq  immediately  to  the 
if  long<=0  %Calculate  msq  of  interest 
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NEARMSQS  (msqind)  =pmmsqs  (pmmsqind)  +floor  (abs  (long/10)  )  ; 
else 

NEARMSQS  (msqind)  =pmmsqs  (pmmsqind)  +35-floor  (long/10)  ; 

end 


end 

end 

end 


Qf  nearby  mSQs*****^*’^***^***-^*^ 

ordflag=0;  %set  'ordered*  flag 
while  ordflag===0 
swapf lag=0; 
for  sortptr= [ 9 : -1 : 2]  ; 

if  NEARMSQS (sortptr)<NEARMSQS (sortptr-1) 
tempmsq=NEARMSQS (sortptr-1) ; 

NEARMSQS (sortptr-1) -NEARMSQS (sortptr) ; 

NEARMSQS ( sortptr ) -tempmsq; 
swapf lag=l ; 

end 

end 

if  swapflag— 0 
ordflag-1; 

end 

end 


all  radiosonde  stations  within  the  9  nearest  msQs’^"^^'^^ 
msqlist-load ( 'msqlist  * ) ; 
listind=0; 
msq-0 ; 

NEARSTATS=[] ; 

for  msqind- [ 1 : 9] ;  %go  through  all  near  msqs 

if  NEARMSQS (msqind) ~-0  %skip  the  out-of-range  msqs 

endlist=0 ; 

while  msq'— NEARMSQS  (msqind)  &endlist— 0 

listind=listind+l;  %flip  through  list  of  msqs  until  you 

msq-msqlist (listind, 1 ) ;  %find  the  one  you  want 
if  listind— length (msqlist) &msq~=NEARMSQS (msqind)  %If  we've 
endlist— 1;  %reached  the  end  of  the  list  w/out  finding  our 
listind=0;  %msq,  we  have  to  end  the  while  loop  and  reset 
end  %listind 

end 

if  endlist==0  %only  add  to  NEARSTATS  if  we  found  our  msq 

NEARSTATS=cat ( 1^  NEARSTATS, transpose ( [msqlist (listind, 2 )  :msqlist (listi 
nd, 2 ) +msqlist (listind, 3) -1] ) )  ; 

end  %  add  the  record  numbers  indicated  in  msqlist  to  the 

end  %  list  of  nearby  stations 

end 

if  is empty (NEARSTATS) 

ISOLATED-1; 

end 
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function  CR4calcdist 
global  RLAT  RLONG  NEARSTATS 

%CLIMAREF  Load  climate  data  module,  submodule  2  (CR4calcdist .m) 

%This  submodule  calculates  the  distances  from  the  radar  to  all  the 
%stations  listed  in  NEARSTATS.  It  then  sorts  the  stations  from  nearest 
%to  furthest. 

% INPUT  VARIABLES: 

%RLAT  -  radar  latitude 
%RLONG  -  radar  longitude 

%NEARSTATS  -  list  of  indices  of  stations  in  the  nine  MSQs  nearest 
%  to  the  radar 

%INTERNAL  VARIABLES: 

%a,b,c  -  earth  radius  related  parameters  (see  code  below) 

%az  -  azimuth  angle  from  radar  to  station  (radians) 

%complat  -  90~latitude  for  a  given  station  (deg) 

%comprlat  -  90-latitude  for  the  radar  site  (deg) 

%deg2rad  -  degrees  to  radians  conversion  factor 

%dellat  -  diff  in  lat  between  a  station  and  the  radar  (deg) 

%dellong  -  diff  in  long  between  a  station  and  the  radar  (deg) 

%dists  -  vector  of  distances  between  radar  and  stations  in  NEARSTATS 

%elevs  -  binary  elevation  values 

%fid  -  file  ID  -  ID#  given  to  radio. dbf 

%latraw  -  binary  latitude  value 

%lats  -  vector  of  latitudes  for  all  stations  in  NEARSTATS  (deg) 

%longraw  -  binary  longitude  value 

%longs  -  vector  of  longitudes  for  all  stations  in  NEARSTATS  (deg) 
%namesandels  -  matrix  binary  codes  for  letters  in  all  the  station  names 
%ordflag  -  flag  indicating  near  stations  in  order  (nearest  to  furthest) 
%psi  -  angle  between  radar  and  a  station  radians) 

%rc  -  local  radius  of  curvature 

%sortptr  -  pointer  used  to  track  stations  in  bubble  sort 
%stat  -  index  used  to  point  to  stations  within  NEARSTATS 
%status  -  dummy  variable  for  I/O 

%swapflag  -  flag  used  in  sort  to  indicate  a  swap  of  positions  has 
happened 

%tempstat  -  temporary  holder  used  in  sort 

%toplO  -  The  number  of  stations  to  keep  (either  10,  or  if  there  aren’t 
%  that  many,  the  number  that  there  are 

%xr,yr,zr  -  3-D  cartesian  coords  of  the  radar (neglecting  radar  altitude) 
%xs,ys,zs  -  3-D  cartesian  coords  of  station  (neglecting  station 
altitude ) 


%0UTPUT  VARIABLES: 

%NEARSTATS  -  see  INPUT  VARIABLES  above,  BUT  NOW  it: 

%  1 .  is  sorted  nearest  to  furthest 

%  2.  contains  distances  from  radar  to  station  in  the 

%  2nd  column 

%  3.  contains  the  ascii  values  of  the  characters  in  the 

%  station  names  in  the  remaining  29  columns 

%  4 .  by  the  time  it’s  output,  it  contains  only  the 

%  10  nearest  stations 


deg2rad=2^pi/360 ;  %degrees  to  radians  conversion  factor 
fid=f open (’ radio . dbf ',’ r ’) ;  %open  file  for  read-only 
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for  stat- [1: length (NEARSTATS) ]  ; 

%First  read  the  lat/longs  for  each  station  in  NEARSTATS 
status-fseek(fid,  642+ (NEARSTATS ( stat ) -1 ) * 94 4  +  38 ,  *bof ' ) ;  %positions 

%pointer  to  start  of  latitude  in  each  record 
latraw=transpose ( fread ( f id, 6) ) ;  %read  raw  (binary)  latitude 
lats (stat ) =str2num (char (latraw) ) ;  %convert  to  number  and  store 

status=fseek(fid,  642+ (NEARSTATS ( stat ) -1 ) * 94 4  +  4 4 ,  ’bof ’ ) ;  %Do  same 

longraw=transpose ( fread ( f id, 7 ) ) ;  %for  longitude 

longs (stat) =-str2num (char (longraw) ) ;  %NOTE:  Data  uses  negative  long. 

%values  to  denote  East  Long.,  but  this  model  uses  negative 

%to  denote  West  Long,  to  be  consistent  with  the  National  Climatic 

%Data  Center  (NOAA)  (hence  the  negative  sign) 


%Next,  calculate  the  local  radius  of  curvature,  RC,  in  the  direction 
%from  the  radar  to  the  station  in  question 
if  sign (RLAT) ==sign (lats (stat ) ) 

dellat^lats (stat ) “RLAT;  %Calc  diff  between  lat/long  of  station  and 
else 

dellat=sign (lats (stat ))* (abs (lats (stat )) tabs (RLAT) ) ;  %if  the  two 
%points  straddle  the  equator,  we  just  add  the  two  abs  val  lats 

end 

if  sign (RLONG) ==sign (longs ( stat ) )  %check  for  longs  straddling  either 

%the  prime  meridian  or  the  int ’ 1  date  line 
dellong=longs (stat ) -RLONG;  %radar  site 
else 

if  abs (RLONG) <90  %i.e.  close  to  zero 

dellong=sign (longs (stat) ) * (abs (RLONG) tabs (longs (stat) ) ) ; %sites 

%straddling  zero  deg 

else 

dellong=sign (longs (stat) ) * (360-abs (RLONG) -abs (longs (stat) ) ) ; 

%sites  straddling  180 

end 

end 

az=pi/2-atan2 (dellat , dellong) ;  %Calc  az  from  radar  to  station 
a=6378.139;  %max  earth  radius  (km) 
b=6356.750;  %min  earth  radius  (km) 
c=a'^2/b^2-l;  %intermediate  value 

rc=a'^2/  (b^sqrt  (l+c^cos  (RLAT*deg2rad)  ^^2)  *  (l+c'^cos  (RLAT*deg2rad)  ^^2^003  (az) 

"2)  )  ; 

%""-radius  of  curvature  (to  be  used  in  place  of  earth’s  radius) 

%Next,  calculate  the  surface  distance  from  the  radar  to  the  station 
complat=90-lats (stat ) ;  %compliment  of  lat(for  spherical  coord  system) 
xs=rc*sin  (complat’^deg2rad)  ^cos  (longs  ( stat )  ’^deg2rad)  ;  %  convert  to 
cartesian  coords 

ys=rc*sin (complat*deg2rad) *sin (longs (stat) *deg2rad)  ; 
zs=rc^cos (complat *deg2rad) ; 

comprlat=90-RLAT;  %compliment  of  lat(for  spherical  coord  system) 
xr=rc*sin  ( comprlat’^deg2rad)  ^cos  (RL0NG*deg2rad)  ;  %convert  to  cartesian 
coords 

yr=rc^sin (comprlat*deg2rad) *sin (RL0NG*deg2rad) ; 
zr=rc*cos  (comprlat'^deg2rad)  ; 
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psi==abs  (acos  {  (xs*xr+ys*yr+zs*zr )  /rc^2)  )  ; 
dists  (stat, 1) =psi*rc; 

end 


%calc  angle  between  radar 
%and  station 


%Concatenate  distances  to  NEARSTATS 
NEARSTATS^cat ( 2 , NEARSTATS , dists ) ; 

%Clear  the  internal  variables 

clear  deg2rad  stat  latraw  lats  longraw  longs  dellat 

clear  dellong  az  a  b  c  rc  complat  xs  ys  zs  comprlat  xr  yr  zr  psi  dists 

%Sort  by  distance  (closest  to  furthest) 
ordflag=0;  %set  ’ordered*  flag 
while  ordflag==0 
swapf lag=0; 

for  sortptr=[size (NEARSTATS, 1 ) :-l:2] ; 

if  NEARSTATS ( sortptr , 2 ) <NEARSTATS ( sortptr-1 , 2 ) 
tempstat=NEARSTATS (sortptr-1,  : )  ; 

NEARSTATS (sortptr-1, : ) =NEARSTATS ( sortptr , : ) ; 

NEARSTATS (sortptr, : ) =tempstat ; 
swapf lag=l ; 

end 

end 

if  swapflag==0 
ordflag=l; 

end 

end 


%Take  the  top  10  nearest  stations 
if  size (NEARSTATS,!) <10 

topl0=size (NEARSTATS, 1) ; 
else 

topl0=10 ; 

end 

NEARSTATS=NEARSTATS  (l:topl0,  :  )  ; 

%Get  the  station  names  and  elevations 
for  stat= [1: size (NEARSTATS, 1) ] ; 

status=f seek (fid,  642+ (NEARSTATS (stat) -1) *944  +  8,  ’bof ’ ) ;  %Get 
name=transpose (fread (fid, 29) ) ;  %station  names  as  well 

status=f seek (fid, 642+ (NEARSTATS (stat ) -1) *944+51, ’bof ’ ) ; 
elevraw=transpose ( f read ( fid, 4 ) ) ; 

elev=double (num2str (round ( str2num (char (elevraw) ) *3 . 28084 ) ) ) ; 
if  length (elev) <5 

for  zros= [1 : 5-length (elev) ] 
elev=[32  elev] ; 

end 

end 

namesandels ( stat, : ) = [name  elev] ; 

end 

NEARSTATS=cat ( 2 , NEARSTATS , namesandels )  ; 

%close  file,  radio. dbf 
status=f close (fid)  ; 
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function  CRSpickstat 
global  USESTAT 
global  NEARSTATS 

%CLIMAREF  Pick  Station  data  module  (CRSpickstat .m) 

%This  module  allows  the  user  to  choose  a  radiosonde  station  from 
%a  list  of  the  nearest  stations.  Average  N  data  from  this  station 
%will  be  used  to  calculate  the  refractivity  profile 

% INPUT  VARIABLES: 

%NEARSTATS  -  Matrix  consisting  of: 

%  Col  1  -  indices  (in  database)  of  all  stations  in  the  nine  Marsden 
%  squares  surrounding  the  radar 

%  Col  2  -  surface  distances  of  the  stations  from  the  radar 
%  Col  3  -  names  of  the  stations  in  binary  form  (i.e.  each  character  in 
%  the  station  name  is  represented  by  its  (decimal)  ascii  code 

%INTERNAL  VARIABLES: 

%statind  -  index  for  for-next  listing  all  the  stations 

%s  -  string  of  near  station  data  and  tabs  for  formatting  display 

%OUTPUT  VARIABLES: 

%USESTAT  -  vector  containing  the  index (ices)  of  the  station (s)  to  use 
%for  the  N-profile.  May  contain  up  to  three  indices. 

format  short 
disp ( ’  ’ ) 

disp(' Index  Name  Height  (ft  MSL)  Distance 

(km)  M 

s=sprintf ( [ ’ 999 ’  ' \t *  ’Standard  Atmosphere  -  ’  ’ \t ’  ’\t’  ’ 

— ']  )  ; 

disp ( s ) 

for  statind=[l :size (NEARSTATS,  1)  ]  ; 

s=sprintf ( [numSstr (NEARSTATS (statind, 1 ) )  ’ \t ’ 

char (NEARSTATS (statind, 3: 36) )  ’ \t '  ’\t’  numPstr (NEARSTATS ( statind, 2 ) ) ] ) ; 

disp(s) ; 

end 

disp  (  ’  ’  ) 

disp ( ’Average  data  from  one, two, or  three  of  these  radiosonde  stations 
must  ’  ) 

disp('be  used  to  construct  the  refractivity  profile.  They  are  listed 
in  ’  ) 

disp (’from  nearest  to  furthest.  You  must  select  the  station (s)  you 
wish ’ ) 

disp (’to  use.  Be  sure  to  consider  topography,  and  not  only  distance 
when ’ ) 

disp (’making  your  decision  (e.g.  may  want  to  choose  a  station  halfway 
between’ ) 

disp(’radar  and  target).  Note:  If  you  elect  to  use  STANDARD 
ATMOSPHERE, ’ ) 

disp (’you  may  not  select  any  other  stations.’) 

USESTAT=[]  ; 

while  isempty (USESTAT) 
disp  (  ’  ’ ) 

disp (’Enter  the  indexes  of  up  to  THREE  stations  (if  you  enter  more 
than’ ) 

USESTAT=input ( ' one,  enclose  the  numbers  in  square  brackets,  e.g  [123 
321]:  ’); 
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if  isempty (USESTAT) | length (USESTAT) >3 

disp(* Enter  1,2,  or  3  station  index  numbers') 

USESTAT=[] ; 

end 

if  length (USESTAT) >1 
saplus=0 ; 

for  i= [ 1 : length (USESTAT) ]  %Check  to  make  sure  the  STD  ATM 
if  USESTAT (i) ==999  %selection  wasn't  chosen  along  with 
saplus=l;  %real  stations 

end 

end 

if  saplus  %If  STD  ATM  was  chosen  along  with  other  stations... 
disp  (  '  ' ) 

disp('You  may  not  use  STANDARD  ATMOSPHERE  along  with  other 
station  data') 

disp  ('It  must  be  used  by  itself) 

USESTAT=[] ; 

end 

if 

( length  ( USESTAT )==2  5cUSESTAT(l)==USESTAT( 2)  )  |  ( length  (USESTAT)  ==3&  (USESTAT 
( 1 ) ==USESTAT (2)1 USESTAT ( 1 ) ==USESTAT ( 3 ) | USESTAT ( 2 ) ==USESTAT ( 3 ) ) ) 

%Check  to  make  sure  duplicate  number  wasn't  entered 
disp  (  '  ' ) 

disp ('You  entered  duplicate  numbers') 

USESTAT=[] ; 

end 

end 

end 
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function  CRGloadNdata 
global  NS  NIK  SEL 
global  USESTAT  MON 

%CLIMAREF  Load  climate  data  module,  submodule  3  (CRGloadNdata .m) 

%This  submodule  loads  Ns  and  Nik  from  the  database  for  the  station  of 
%interest 

% INPUT  VARIABLES: 

%USESTAT  -  index  of  station  of  interest 
%MON  -  month  of  interest 

%INTERNAL  VARIABLES: 

%delmlk  -  m-gradient  (numeric) 

%delmlkraw  -  raw  (binary)  m-gradient 
%delnlk  -  n-gradient  (ngradient=mgradient-156) 

%fid  -  file  id 

%monstart  -  position  in  record  of  the  start  of  the  data  for  the  month  of 
interest 

%nlk  -  scalar  value  of  lOR  at  1km  AGL  for  one  particular  station 

%ns  -  scalar  value  of  ns  for  one  particular  station  -  is  cat ' d  onto  the 

NS  vector  eventually 

%nsraw  -  raw  (binary)  surface  refractivity 

%sel  -  scalar  value  of  the  station  elevation  for  one  particular  station 
%selraw  -  raw  (binary)  surface  elevation 
%stat  -  index  to  USESTAT  in  for/next  loop 
%status  -  dummy  variable  for  i/o  operations 


%OUTPUT  VARIABLES: 

%NS  -  surface  index (ices)  of  refraction 

%N1K  -  index (ices)  of  refraction  at  1km  above  surface 

%SEL  -  station  elevations  (km) 


if  USESTAT-=999 

fid=f open (’ radio . dbf r *) ;  %open  file  for  read-only 

NS=[]  ; 

N1K=[]  ; 

SEL=[]  ; 

for  stat=[l:length{USESTAT) ] 

status=f seek ( fid,  642+ (USESTAT (stat) -1 ) *94  4  +  51,  * bof ’ ) ;  %positions 

%pointer  to  start  of  elevation  in  each  record 
selraw=transpose ( f read ( f id, 6) ) ;  %read  raw  (binary)  elevation 
sel= . 001*str2num (char ( selraw) ) ;  %convert  to  number  (km)  and  store 
SEL=cat (2, SEL, sel) ;  %tack  onto  SEL  vector 


monstart-642+ (USESTAT (stat) -1) *944+55+ (MON-1) *74;  %pointer  to 
start  of  data 


%for  each  month  in  the  record 


status=f seek ( fid, monstart+12 ,' bof ') ;  %pointer  to  start  of  Ns  in 

%each  record 

nsraw=transpose (fread(fid, 3) ) ;  %read  raw  (binary)  Ns 
ns=str2num { char (nsraw) ) ;  %convert  to  number  and  store 
NS=cat ( 2 , NS , ns ) ; 

status=fseek (fid,monstart+15, 'bof ' ) ;  %pointer  to  start  of  Ns  in 

%each  record 

delmlkraw=transpose ( f read ( f id, 3 ) ) ;  %read  raw  (binary)  Ns 
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delnilk=str2num  (char  (delmlkraw)  )  ;  %convert  to  number  and  store 
delnlk=delmlk-156;  % 

nlk=ns+delnlk;  %Calculate  N  at  Ik  from  surface  N  and  gradient 
NlK-cat (2,NlK,nlk) ; 

end 

status=fclose (fid) ; 
else 

SEL=0; 

NS=313;  %User  chose  standard  refraction 
N1K=271.0612; 

end 

%disp{[’NS=’  num2str(NS)  '  N1K=*  num2str (NIK) ] )  %for  testing  only 
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function  CR7interp 

global  USESTAT  NEARSTATS  RLAT  RLONG  NS  NIK  SEL 

%%CLIMAREF  Interpolate  module  (CR7interp .m)  ,  ,  ,  o  o 

%This  submodule  interpolates  NS  and  NIK  from  data  recorded  for  2  or  3 

%stations 


%INPUT  VARIABLES: 

%USESTAT  -  indices  of  stations  of  interest 

%NEARSTATS  -  1.  contains  indexes  of  nearest  ten  stations  in  col  1 

%  2 .  contains  distances  from  radar  to  station  in  col  2 

%  3.  contains  the  ascii  values  of  the  characters  in  the 

%  station  names  in  the  remaining  29  columns 

%  4 .  is  sorted  nearest  to  furthest 

Q. 

"O 

%RLAT  -  radar  latitude 

%RLONG  -  radar  longitude  ^  . 

%USESTAT  -  vector  of  station  indexes  to  be  used  in  interpolation 
%NS  -  surface  lOR  vector  (for  both,  or  all  three  stations)  -  indexed 
%  to  match  USESTAT 

%N1K  -  IK  above  surface  lOR  vector  (for  both,  or  all  three  stations) 
%  --  indexed  to  match  USESTAT 

%SEL  -  surface  elevations  of  stations 


%INTERNAL  VARIABLES: 

%a,b,c  -  earth  radius  related  parameters  (see  code  below)  ^ 

%a,b, cns, cnlk, csel  -  direction  numbers  for  the  line  used  in  the  line- 
interp 

%Ans, Anlk, Bns, Bnlk, C  -  direction  numbers  for  plane  used  in  plane-interp 
%az  azimuth  angle  from  radar  to  station  (radians) 

%ce  -  intermediate  constant  for  N-profile  calculation  (in  this  case  for 
normalization  of  NS  and  NIK 

%complata  ”  90“latitude  for  first  station  site  (deg) 

%complatb  -  90-latitude  for  second  station  (deg) 

%d23a,d23b  -  two  possible  alternative  values  for  the  dist  from  station  2 
to  station  3  (only  one  is  correct) 

%dellat  -  diff  in  lat  between  a  station  and  the  radar  (deg) 

%dellong  -  diff  in  long  between  a  station  and  the  radar ^  (deg) 

%dist  -  index  in  distance-finding  for/next  loop  to  station  a 

%distb  -  index  in  distance-finding  for/next  loop  to  station  b 

%dists  -  vector  of  distances  between  stations  (1-2)  or  ( 1-2 , 2-3, 3-1 ) 

%Dns,Dnlk  -  parameter  which  locates  plane  in  space 

%latraw,  longraw  -  binary  version  of  lat  &  long  for  station 

%lats,  longs  -  lats  and  longs  for  stations  in  USESTAT 

%nlkmsl  -  Nik  normalized  to  sea  level  (i.e.  N  at  1km  MSL) 

%nsi,nlki  -  interpolated  values  of  ns  and  nlk  still  normalized  to  MSL 
%nsmsl  -  Ns  normalized  to  sea  level  (i.e.  N  at  0km  MSL) 

%numdists  —  number  of  distances  between  stations  to  be  calculated 
%psi  -  angle  between  station  a  and  station  b  radians) 

%rc  -  local  radius  of  curvature 

%rdists  -  vector  of  distances  from  radar  to  each  station  (r-l,r-2)  or 
(r-1, r-2, r-3) 

%stat  -  index  to  USESTAT  in  for/next  loop 
%status  -  dummy  variable  for  i/o 

%xl, yl, x2, y2, x3, y3  -  station  coordinates  (for  interpolation  purposes) 
%xa,ya,za  -  3-D  cartesian  coords  of  station  a  (neglecting  radar 
altitude) 

%xb,yb,zb  -  3-D  cartesian  coords  of  station  b (neglecting  station 
altitude ) 
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%xp  -  X'-coordinate  of  point  on  projection  of  line  (onto  2d  plane)  that’s 
closest  to  the  radar 

%xr,yr  -  radar  coordinates  (for  interpolation  purposes) 

%y3a^y3b  -  two  possible  alternative  values  for  y3  (only  one  is  correct) 


%OUTPUT  VARIABLES: 

%NS  -  interpolated  index  of  refraction  at  interpolated  surface 
%N1K  -  interpolated  index  of  refraction  at  1km  above  interpolated 
%  surface  (MSL) 

%SEL  -  interpolated  station  elevation  (km)  (interpolated  surface) 


deg2rad=2^pi/360 ;  %degrees  to  radians  conversion  factor 

fid=f open (’ radio . dbf r ’) ;  %open  file  for  read-only 

%%%%%%Read  the  lat/longs  for  each  station  in  NEARSTATS%%%%% 
for  stat- [1: length (USESTAT) ] ; 

status=fseek(fid, 642+ (USESTAT (stat) -1) *944+38, 'bof ’ ) ;  %positions 

%pointer  to  start  of  latitude  in  each  record 
latraw=transpose ( fread (f id, 6) ) ;  %read  raw  (binary)  latitude 
lats ( stat ) =str2num (char (latraw) ) ;  %convert  to  number  and  store 

status=fseek(fid, 642+ (USESTAT (stat ) -1) *944+44, 'bof ’ ) ;  %Do  same 

longraw=transpose (fread (fid, 6) ) ;  %for  longitude 

longs (stat) =-str2num (char (longraw) ) ;  %NOTE:  Data  uses  negative  long. 

%values  to  denote  East  Long. ,  but  this  model  uses  negative 

%to  denote  West  Long,  to  be  consistent  with  the  National  Climatic 

%Data  Center  (NOAA)  (hence  the  negative  sign) 

end 

status=f close (fid) ; 

%%%%%Calc  distances  between  the  stations%%%%%% 
if  length (USESTAT)=-2 

numdists=l;  %If  there ’ re  two  stations  calc  one  distance 
else 

numdists=3;  %If  there ’re  three  stations  calc  three  distances 

end 

for  dist= [ 1 : numdists] 

distb=dist+l ;  %indexing  the  second  station 
if  distb==4 
distb=l ; 

end 

%calculate  the  local  radius  of  curvature,  RC,  in  the  direction 
%from  the  first  station  to  the  second  station 
if  sign (lats (dist ) ) ==sign (lats (distb) ) 

dellat=lats (distb) -lats (dist ) ;  %Calc  diff  between 
%lat/long  of  station  a  and  station  b 
else 

dellat=sign (lats (distb) ) * (abs (lats (distb) ) tabs (lats (dist ) ) ) ;  %if 
the  two 

%points  straddle  the  equator,  we  just  add  the  two  abs  val  lats 

end 

if  sign  (longs  (dist ))  ===sign  (longs  (distb)  )  %check  for  longs  straddling 
either 
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%the  prime  meridian  or  the  int ’ 1  date  line 
dellong=longs (distb) -longs (dist) ;  %radar  site 
else 

if  abs ( longs (dist )) <90  %i.e.  close  to  zero 

dellong=sign (longs (distb) ) * (abs (longs (dist) ) +abs (longs (distb) ) ) ;% sites 

%straddling  zero  deg 

else 

dellong=sign (longs (distb) ) * (360-abs (longs (dist) ) - 
abs (longs (dist) ) ) ; 

%sites  straddling  180 

end 

end 

az=pi/2-atan2 (dellat , dellong) ;  %Calc  az  from  radar  to  station 
a=6378.139;  %max  earth  radius  (km) 
b=6356.750;  %min  earth  radius  (km) 
c==a^2/b'^2-l ;  %intermediate  value 

rc=a^2/  (b^sqrt  (l  +  c*cos  (lats  (dist)  *deg2rad)  ^^2)  (l  +  c*cos  (lats  (dist)  *deg2ra 
d)  ^2^cos  (az)  ^"2)  )  ; 

%^-radius  of  curvature  (to  be  used  in  place  of  earth’s  radius) 

%Next,  calculate  the  surface  distance  from  the  first  station  (a)  to 
%the  second  station  (b) 

complatb=90-lats (distb) ;  %compliment  of  lat(for  spherical  coord 
system) 

xb=rc*sin  ( complatb*deg2rad)  ’^cos  (longs  (distb)  ’^deg2rad)  ;  %convert  to 
cartesian  coords 

yb^rc^sin  (complatb*deg2rad)  ^sin  (longs  (distb)  ^deg2rad)  ; 
zb=rc*cos (complatb*deg2rad)  ; 

complata=90-lats (dist ) ;  %compliment  of  lat(for  spherical  coord 
system) 

xa^rc^sin  (complata’^deg2rad)  ’^cos  (longs  (dist)  *deg2rad)  ;  %convert  to 
cartesian  coords 

ya=rc*sin ( complata*deg2rad) *sin (longs (dist) *deg2rad) ; 
za=rc^cos (complata*deg2rad) ; 

psi=abs  (acos  (  (xb*xa+yb*ya+zb^za) /rc'^2 )  )  ;  %calc  angle  between  radar 

%and  station 

dists  (dist )  =psi^ res¬ 
end 

%also  extract  distances  from  radar  to  stations  from  NEARSTATS 
for  stat= [ 1 : length (USESTAT) ]  %go  through  list  of  2  or  3  stations 
for  nstat= [ 1 : size (NEARSTATS, 1 ) ]  %go  through  list  of  10  NEARSTATS 
if  USESTAT  (stat)=NEARSTATS  (nstat,  1) 
rdists (stat ) =NEARSTATS (nstat, 2) ; 

end 

end 

end 

%%%%%%%%Construct  triangle  using  three  distances:  (r-1, r-2, 1-2) %%%%%%% 
%%%%%%%%This  routine  is  used  whether  you  have  two  stations  or  three%%% 
%%%%%%%%This  routine  essentialy  flattens  out  the  surface  of  the%%%%%%% 
%%%%%%%%earth  in  order  to  do  a  simple  planar  or  linear  interpolation%% 
xr=0;  %establish  radar  as  center  of  coord  system 
yr=0; 

xl=rdists ( 1 ) ;  Iposition  station  2  on  the  x  axis  to  the  right  of  radar 
yl=0;  %Note:  absolute  orientation  of  this  system  doesn't  matter  for 
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%  interp,  only  relative  orientation 

x2=  (rdists  (2)  ^2-dists  (1)  '^2+rdists  (1)  -^2)  /  (2^rdists  (1)  )  ;  %Calc  x2,  y2 
y2=sqrt  (rdists  (2 ) '^2-x2'^2 )  ;  %Based  on  the  other  two  points  —  in  addition 
%to  its  arbitrary  orientation,  this  triangle  may  also  be  flipped 
% (mirror  image).  This  doesn't  matter  either  for  the  interpolation 

%%%%%%%Locate  third  station  (another  triangle)  if  there  is  one%%%%%%% 
%%%%%%%This  is  used  only  for  the  planar  (3-station)  interpolation%%%% 
if  length (USESTAT) ==3  %locate  third  station  if  there  is  one 
x3=  (rdists  (3)  '^2-dists  (3)  '^2+rdists  (1)^2)/  (2*rdists  (1)  )  ; 
y3a=sqrt (rdists (3) ^2-x3^2) ;  %station  three  is  either  above  or 
y3b=-y3a;%  below  the  x-axis 

d23a=sqrt ( (x2-x3) ^2+ (y2-y3a) ^2 )  ; 
d23b=sqrt { (x2-x3) ^2+ (y2-y3b) ^2)  ; 

if  abs (d23a“dists (2) ) <abs (d23b-dists (2) )  %To  figure  out  which,  we 
y3=y3a;%  simply  calculate  both,  figure  out  the  two  distances  to 
else  Istation  2,  and  see  which  one  matches  the  actual  distance 
y3=y3b;  %the  closest 

end 

end 


%%%%%%Normalize  NS  &  NIK  (all  stations)  to  Oft  MSL%%%%%%%%% 
ce=log(NS. /NIK) ; 

nsmsl=NS . *exp (-ce . * (-SEL) ) ;  %calc  N  at  sea  level 
nlkmsl=NS , *exp (-ce . *  (-SEL+1) ) ;  %calc  N,  1km  above  sea  level 


%%%%%%Linear  Interpolation  (for  2  stations) 

%construct  3  sets  of  lines  in  3d  space  -  one  using  nsmsl  for  z, 
%another  using  nlkmsl  for  z,  and  the  third  using  SEL  for  z 
if  length (USESTAT)==2 

a=x2-xl;  %parametric  coefficients  for  lines 
b=y2-yl; 

cns=nsmsl (2) -nsmsl (1) ; 
cnlk=nlkmsl (2) -nlkmsl (1); 
csel=SEL(2) -SEL{1) ; 


%Find  point  on  projection  (of  line  onto  2d  plane)  closest  to  radar 
xp= (a*xr/b+b*xl/a-yl+yr ) / (b/a+a/b) ;  %We  only  need  x-coord--it ' s  a 
line ! 

%Find  z-value  (ns  or  nlk) 

nsi=cns* (xp-xl) /a+nsmsl (1) ;  %Note;  We're  redefining  NS, NIK, &  SEL  as 
nlki=cnlk^ (xp-xl ) /a+nlkmsl ( 1) ;  %scalar  interpolated  values 
SEL=csel’^  (xp-xl)  /a+SEL  ( 1)  ; 
else  %i.e.  there  are  3  stations 

Ans= (y2-yl) * (nsmsl (3) -nsmsl (1) ) - (nsmsl (2) -nsmsl ( 1 ) ) * ( y3-yl ) ; 

Anlk= (y2-yl) * (nlkmsl (3) -nlkmsl (1) ) - (nlkmsl (2) -nlkmsl (1) ) ^ (y3-yl) ; 
Asel-(y2-yl)^(SEL(3)-SEL(l) ) - ( SEL { 2 ) -SEL ( 1 ) )*(y3-yl) ; 

Bns= (nsmsl (2) -nsmsl (1) ) * (x3-xl) - (x2-xl) * (nsmsl (3) -nsmsl (1) ) ; 

Bnlk= (nlkmsl (2) -nlkmsl (1) ) ^ (x3-xl) - (x2-xl) * (nlkmsl (3) -nlkmsl (1) ) ; 
Bsel- (SEL (2) -SEL ( 1) ) * (x3-xl) - (x2-xl) * (SEL ( 3 ) -SEL ( 1 ) ) ; 

C=(x2-xl) * (y3-yl)- (y2-yl) * (x3-xl) ; 

Dns=-Ans*xl-Bns*yl“C^nsmsl (1) ; 

Dnlk=-Anlk*xl-Bnlk*yl-C^nlkmsl (1) ; 

Dsel=-Asel*xl-Bsel^yl-C^SEL(1) ; 


nsi= (-Ans*xr-Bns*yr-Dns) /C;  %Redefine  NS, NIK, &  SEL  as  scalars  with 
nlki= (-Anlk^xr-Bnlk*yr-Dnlk) /C; %the  interpolated  values 
SEL= (-Asel*xr-Bsel*yr-Dsel) /C; 
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end 


%De-normalize  NS  and  Nik  back  to  the  interpolated  SEL 
ce=log (nsi/nlki ) ; 

NS=nsi*exp (-ce^SEL) ;  %calc  N  at  surface 
NlK=nsi^exp  (“ce"^  (SEL+1)  )  ;  %calc  N,  1km  above  surface 
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function  CRSductstats 

global  USESTAT  MON  DYNT  SEL  MULTSTATFLG 
global  GMSB  OHSB  MTSB  MDSB  MFSB  PSB  GMEL  OREL  MTEL 
global  MDEL  MFEL  PEL  P2EL  PSBEL  NS  NIK 
global  PHs  PMs  PHe  PMe  PHn  PMn 

%CLIMAREF  Load  climate  data  module,  submodule  4  (CRSductstats .m) 

%This  submodule  loads  ducting  statistics  from  the  database  for 
%the  station  of  interest 

% INPUT  VARIABLES: 

%USESTAT  -  index  of  station  of  interest 
%MON  -  month  of  interest 

%DYNT  -  day,  night,  or  both  readings  (l=day,  0=night,  2=both) 

% INTERNAL  VARIABLES: 

%fid  -  file  id 

%monstart  -  position  in  record  of  the  start  of  the  data  for  the  month  of 
interest 

%p2elxl00  -  prob  >1  elev  duct  occurs  (%xl00) 

%pel00  -  %  of  OOOOZ  radiosonde  readings  in  which  elev.  ducts  occurred 

%pell2  -  %  of  1200Z  radiosonde  readings  in  which  elev.  ducts  occurred 

%psb00  -  %  of  OOOOZ  radiosonde  readings  in  which  surf,  ducts  occurred 

%psbl2  -  %  of  1200Z  radiosonde  readings  in  which  surf,  ducts  occurred 

%psbelxl00  -  prob  surf-based  and  elev  duct  both  occur  simultaneously 
(%xl00) 

%status  -  dummy  variable  for  i/o  operations 
%OUTPUT  VARIABLES: 

% (heights  in  kilometers  unless  otherwise  specified) 

%GNEL  -  elevated  duct  N-unit  gradient 
%GNSB  -  surface  based  duct  N-unit  gradient 
%NDEL  -  elevated  duct  N-unit  deficit 
%NDSB  -  surface  based  duct  N-unit  deficit 
%NFEL  -  elevated  duct  trapping  frequency 
%NFSB  -  surface  based  duct  trapping  frequency 
%NTEL  -  elevated  duct  thickness 
%NTSB  -  surface  based  duct  thickness 
%OHEL  -  elevated  duct  optimum  height 
%OHSB  -  surface  based  duct  optimum  height 
%P2EL  -  probability  of  >1  elevated  duct 

%PEL  -  elevated  duct  percent  chance  of  occuring  (day,  night,  or  avg) 

%PHe  -  vector  of  key  heights  (km)  in  elevated  duct  profile 

%PHn  -  vector  of  key  heights  (km)  in  no-duct  profile 

%PHs  -  vector  of  key  heights  (km)  in  surface  duct  profile 

%PMe  -  vector  of  M  values  corresponding  to  heights  in  PHs 

%PMn  -  vector  of  M  values  corresponding  to  heights  in  PHs 

%PMs  -  vector  of  M  values  corresponding  to  heights  in  PHs 

%PSB  -  surface  based  duct  percent  chance  of  occuring  (day,  night  or  avg) 

%PSBEL  -  probability  of  surface-based  AND  elevated  duct  occuring 

%Initialize  variables 
PHs=[] ; 

PMs- [ ] ; 

PHe-[]  ; 

PMe- [ ] ; 

PHn- [ ] ; 

PMn- [ ] ; 

%Get  Ducting  Data 


A-26 


if  USESTAT-=999&~MULTSTATFLG 

fid=fopen( 'radio.dbf , ) ;  %open  file  for  read-only 

monstart=642+ (USESTAT-1) *944+55+ (MON-1) *74/  %pointer  to  start  of 

%data  for  each  month  in  the  record 

status=fseek(fid,monstart+18,  ^bof ' ) ;  %pointer  to  start  of  duct  data 

ductdatraw=transpose (fread(fid, 56) ) ;  %read  entire  list  of  raw 

% (binary)  duct  statistics 

%Now,  break  it  out  and  convert  the  binary  data  to  actual  numbers 
GMSB=str2num(char (ductdatraw(l:4) ) ) ;  %convert  gmsb  from  binary,  to 
%  a  string,  then  to  a  number 

0HSB=str2num(char (ductdatraw(5 : 7) ) ) /lOOO;  %  x/lk  to  convert  to  km 
MTSB=str2num(char (ductdatraw ( 8 :10)))/1000; 

MDSB=str2num (char (ductdatraw (11:13) ) )  ; 

MFSB=str2num (char (ductdatraw ( 14 : 17)  )  )  ; 
psbl2=str2num (char (ductdatraw (18: 20) )  )  ; 
psb00=str2num (char (ductdatraw (21:23))); 

if  DYNT=0  %determine  which  parameter  to  use  depending  on  user 
PSB=psb00;  %input  of  day,  night,  or  both 
elseif  DYNT==1 
PSB=psbl2; 
else 

PSB= (psb00+psbl2) /2; 

end 

GMEL=str2num (char (ductdatraw (24:27))); 

OHEL=str2num (char (ductdatraw (28:31) ) ) /lOOO; 

MTEL=str2num (char (ductdatraw (32:34)  ))/1000; 

MDEL=str2num (char (ductdatraw (35 : 37 ) ) )  ; 

MFEL=str2num (char (ductdatraw (38 :41))); 
pel 12=str2num (char (ductdatraw (42:44))); 
pel00=str2num (char (ductdatraw ( 45 : 47 )  )  )  ; 

if  DYNT==0  %determine  which  parameter  to  use  depending  on  user 
PEL==pel00;  %input  of  day,  night  or  both 
elseif  DYNT==1 
PEL-pell2; 
else 

PEL=(pel00+pell2) /2; 

end 

p2elxl00=str2num (char (ductdatraw (48:51))); 
psbelxl00=str2num (char (ductdatraw (52:56)  )  )  ; 

P2EL-p2elxl00/100; 

PSBEL=psbelxl00/100; 

status=f close (fid) ; 

%Calculate  duct  heights  and  associated  M  values 
MlK-NlK+1* (le6/6378) ; 

%  Surface  Duct  Calculations 
PHs(l)=0;  %all  heights  in  km 
PMs (1)=NS; 

PHs (2)=OHSB; 

PMs (2) =NS+GMSB*OHSB; 

PHs (3) =MTSB; 

PMs (3)=PMs (2)-MDSB; 

PHs (4)=PHs (3)+l/ 

PMs (4)=PMs (3)+ (MIK-NS) ; 
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%  Elevated  Duct  Calculations 
PHe(l)=0;  %all  heights  in  km 
PMe (1)=NS; 

PHe (2)=OHEL; 

PMe ( 2 ) =NS+GMEL*OHEL ; 

PHe (3 ) =OHEL- (MDEL/GMEL) +MTEL; 

PMe(3)=PMe(2)-MDEL; 

PHe(4)=PHe{3)+l; 

PMe(4)=PMe(3)+(MlK-NS) ; 

%  No-Ducting  Calculations 
PHn ( 1 ) =0 ; 

PMn ( 1 ) =NS ; 

PHn(2)=l; 

PMn(2)=MlK; 

if  ~PHs(2)  %Happens  if  the  prob  of  surf  duct  is  zilch,  so  we 
PHs (2) =PHs (4) ;  %just  calculate  the  two  levels  needed  for 
PMs (2 ) =PMs (4 ) ;  %the  exponential  model 
PHs(3:4)=[] ; 

PMs ( 3 : 4 ) = [ ] ; 

end 

if  ~PHe(2)  %Happens  if  the  prob  of  elev  duct  is  zilch,  so  we 
PHe (2) =PHe (4 ) ;  %just  calculate  the  two  levels  needed  for 
PMe (2) =PMe (4 ) ;  %the  exponential  model 
PHe{3:4)=[] ; 

PMe ( 3 : 4 ) = [ ] ; 

end 

else  %Standard  atmosphere  or  interpolation  -  no  ducting 
GMSB=0; 

OHSB=0; 

MTSB=0; 

MDSB=0; 

MFSB=0; 

PSB=0 ; 

GMEL=0; 

OHEL=0; 

MTEL=0; 

MDEL=0; 

MFEL=0; 

PEL=0 ; 

P2EL=0; 

PSBEL=0; 

%  No  Ducting  Profile  Calculations 
M1K=N1K+1* (le6/6378) ; 

PHn ( 1 ) =0 ; 

PMn ( 1 ) =NS ; 

PHn ( 2 ) =1 ; 

PMn(2)=MlK; 

end 
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function  CR9raytrace 

global  DUCTFLG  PHs  PHe  PHn  PMs  PMe  PMn  RAZ  REL  RLAT  SEL  THT  TRG 
global  ALFO  ALFG  BETAES  RO  HEIGHTS  THTAPP  THTERR  THTERRIOO  TRGAPP 
global  TRGERR  TRGERRIOO  DUCTHTS  ERRMSG 

%CLIMAREF  Raytrace  with  ducting  module  (CR9raytrace_d.m) 

%This  module  uses  the  refractivity  data  from  the  selected  radiosonde 
%station  to  trace  the  bending  of  a  ray  from  the  radar  to  the  given 
%target  whether  or  not  a  duct  exists 

% INPUT  VARIABLES: 

%DUCTFLG  -  flag  indicating  what  type^  if  any,  ducting  will  be  taken  into 
account 

%PHs,PHe,PHn  -  M-profile  heights 

% PMs, PMe, PMn  -  M-profile  M-values  into  account  ( 0=none, l=sf c, 2=elev) 

%RAZ  -  radar  azimuth  (direction  it's  pointing) 

%REL  -  radar  elevation  (AGL) 

%RLAT  -  radar  latitude 

%SEL  -  radiosonde  station  height  (MSL) 

%THT  -  target  height  (AGL)  (actual)  (km) 

%TRG  -  target  range  (actual)  (km) 

%INTERNAL  VARIABLES: 

%a,b,c  -  intermediate  variables  for  radius  of  curvature  calc 
%A,B,C,delz  -  intermediate  values  in  calculating  height  of  extrema 
%alfga,alfgb  -  intermediate  variables  for  calculating  ALFG 
%a0  -  estimate  of  initial  take-off  angle  required  to  hit  target 
%aOhiB  -  highest  aO  of  a  trace  that  reached  target  ht  but  missed  tgt 
%a01astgood  -  last  aO  estimate  that  was  "in  the  clear"  —  used  for 
%  splitting  the  difference 

%a01oB  -  lowest  aO  of  a  trace  that  reached  target  ht  but  missed  tgt 
%aOmagtoolo  -  lower  limit  of  takeoff  angle  when  radar  and  target  are  in 
%  the  duct  together.  Takeoff  angle  <=  this  will  result  in 

%  trace  never  reaching  the  target  height 

%aOmagtoohi  -  upper  limit  of  takeoff  angle  when  radar  and  target  are  in 
the 

%  duct  together 

%aOtoohi  -  an  aO  estimate  above  aOtoohi  will  either  put  the  trace  in  a 
%  duct  unnecessarily  or  cause  the  ray  to  overshoot  the  target 

%  on  a  downward  slope 

%aOtoolo  -  an  aO  estimate  below  aOtoolo  will  either  put  the  trace  into 
%  a  duct  unnecessarily  or  cause  the  ray  to  "hit  dirt" 

%  unnecessarily 

%angtol  -  tolerance  for  subtense  matching  (constant=le-6) 

%bend  -  index  into  profile  heights,  ph 
%betal  -  subtense  over  current  layer 

%beta43, rc43, trg43  -  subtense, radius  of  curvature, and  target  range 
%  modified  for  4/3  earth 

%betae  -  total  subtense  up  to  this  point  in  the  estimated  trace 
%betat  -  actual  subtense  between  radar  and  target 
%botductht  -  height  of  bottom  of  duct 

%bothinduct  -  flag  indicating  both  radar  and  target  are  in  the  duct 
%botrng, toprng  -  defines  levels  at  top  and  bottom  of  a  profile  layer 
%cantfind  -  flag  indicating  the  estimation  is  not  working  to  find  the 
%  target,  probably  because  of  ducting  effects 

%ce  -  intermediate  constant  in  calculation  of  N 

%clostsub  -  when  trace  has  reached  target  height  multiple  times, 

%  dost  sub  is  the  index  into  hitBETAES ,  and  hitdBdaO, 
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%  indicating  which  is  the  closest  to  the  target 

%cosial, cosiaO, cosia2  -  cosine  of  al,a0,a2 

%dbldaO  -  incremental  subtense  change  with  initial  el  angle 
%dbdaO  -  total  subtense  change  with  initial  el  angle 
%deg2rad  -  degres  to  radians  conversion 
%dh  -  height  increment 

%gam  -  intermediate  var  --  n-gradient  in  layer 
%h  -  height  vector,  contains  heights  of  all  defined  levels 
%hitbetae  -  subtenses  of  points  at  which  target  height  was  reached 
%hitdBdaO  “  total  change  in  subtense  with  change  in  aO,  stored  for 
%  every  point  at  which  target  height  was  reached 

%hitdirtflg  -  flag  indicating  ray  has  hit  the  ground 

%hithole  -  flag  indicating  target  appears  to  be  in  a  radio  hole  (val=l) 

%  —  duct  steers  all  radio  energy  away  from  it  --  OR  over 

%  the  radio  horizon  (val=2) 

%hittarg  -  flag  indicating  target  hit  (value  indicates  which  hit, 

%  if  more  than  one,  was  the  target) 

%maxalt  -  maximum  altitude  for  defining  refractivity  profile 
%maxminflg  -  indicates  an  extrema  (ray  is  tangent  to  a  shell 
%  concentric  with  the  earth)  has  been  reached 

%Mgrad  ^  M  gradient,  used  for  deteriming  bottom  of  elevated  duct 
%missby  -  subtense  angle  (pos  or  neg)  by  which  the  trace  missed  target 
%N  -  refractivity  at  every  level  from  surface  to  maxalt 
%n0  -  index  of  refraction  at  surface 

%Nl,hl  -  N  and  h  values  for  lower  level  defining  current  layer 
%N2,h2  -  N  and  h  values  for  upper  level  defining  current  layer 
%Nm, hm  -  N  and  h  values  at  extrema 
%nummms  -  extrema  counter 

%outoduct  -  flag  indicating  trace  has  left  the  duct 
%ph,pM  ~  M-profile  heights  and  M  values 
%pN  -  N-profile  N-values 

%psil  -  amount  of  bending  ray  underwent  in  current  layer 

%R0  ”  (rc+SEL)  radius  from  earth  center  to  surface  (station)  elevation 

%rta,rtb  -  intermediate  variables  for  calculating  ALFG 

%stepstohit  -  number  of  steps  it  took  to  get  to  each  point  at  which 

%  target  height  was  reached 

%stoptrace  -  flag  indicating  trace  must  be  stopped  -  value  indicates 
%  reason  for  stop  (see  notes) 

%trgappl  -  range  from  point  1  to  point  2  for  a  single  layer  — 

%  accumulated  to  calculate  TRGAPP,  the  total  apparent  range 

%OUTPUT  VARIABLES: 

%ALF0  -  initial  angle  of  elevation  for  iteration 

%ALFG  -  geometric  elevation  angle  (el  angle  of  straight  line  between 
%  radar  and  target 

%BETAES  -  stored  betaes  for  plotting  later 

%HEIGHTS  -  stored  step  heights  for  plotting  later 

%R0  -  local  earth  radius  of  curvature  +  surface  elevation 

%THTAPP  -  apparent  target  height 

%THTERR  -  height  error 

%THTERR100  -  height  percent  error 

%TRGAPP  -  total  apparent  range 

%TRGERR  -  range  error 

%TRGERR100  ■“  height  percent  error 
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initial  variables  and  values*^***^*^*****^*^^****^ 


inaxnumtraces=30 ; 

deg2rad=2^pi/360;  %degrees  to  radians  conversion  factor 
angtol=le-6;  %tolerance  used  to  compare  angles  in  iteration  (radians) 
dh=.001;  %height  resolution  in  km  (recommend  dh=0. 00001  for  best  results 
maxalt=5+max ( [THT  REL] ) ;  %we  won’t  worry  about  levels  any 
%  higher  than  5km  above  highest  point  —  that'll  cover  ducting 
maxalt=dh*round(maxalt/dh) ;  %round  maxalt  to  nearest  level 

REL=dh'*^round  (REL/dh)  ;  %round  radar  elevation  to  nearest  level 

THT=dh*round(THT/dh) ;  %round  target  height  to  nearest  level 

^Calculate  local  earth  radius  of  curvature 
a=6378.139;  %km 
b-6356.750;  %km 
c=a'"2/b'^2~l; 

rc-a'^2/  (b^sqrt  (l  +  c*cos  (RLAT*deg2rad)  ^"2)  *  (l+c*cos  (RLAT*deg2rad)  ^2^ 
cos  (RAZ’'deg2rad)  ^2)  )  ; 

%define  radius  to  surface 
R0=rc+SEL; 

%Calculate  N-Profile 

switch  DUCTFLG  %determine  which  ducting  profile,  if  any,  to  use 
case  0 

ph=PHn; 
pM^PMn; 
case  1 

ph=PHs ; 
pM=PMs ; 
case  2 

ph=PHe; 

pM==PMe; 

end  %end  of  determine  ducting  profile 

ph=dh*round(ph/dh) ;  %round  profile  heights  to  nearest 
%  multiple  of  dh 

pN=pM-ph* (le6/6378) ; 

pgd= (pN (2 : length (pN) ) -pN (1 : length (pN) -1) ) . / . . . 

(ph (2 : length (ph) ) -ph (1 : length (ph) -1)); 
h=[0:dh:maxalt] ;  %height  vector  all  heights  (hence  levels)  are 
%  defined  w.r.t.  surface  elevation 

for  bend= [1 : length (ph) -2] 

botrng=round(ph(bend) /dh+1) ;  %calculates  bottom  and  top  levels  of 
toprng=round(ph(bend+l) /dh+1) ;  %layer  —  'round'  ensures  integer 
N (botrng: toprng) =pN (bend) + (h (botrng: toprng) -h (botrng) ) *pgd (bend) ; 
end  %end  of  loop  through  all  "bends" 

ce=log (pN (length (pN) -1) /pN (length (pN) ) ) ;  %intermediate  constant 
botrng=round (ph (length (ph) -1) /dh+1) ;  %'round'  simply  ensures  answer 
toprng=round(raaxalt/dh+l) ;  %is  integer  type — no  rounding  really  tajces 
place 

N (botrng: toprng) =pN (length (pN) -1) * . . . 

exp (-ce* (h (botrng : toprng) -h (botrng) ) ) ; 

nO=l+N(l) *le-6;  lindex  of  refraction  at  surface 
%Calculate  target  subtense 

betat=acos  {  {  (RO+THT)  ^2+  (RO+REL)  ■'2-TRG''2)  /  (2*  (RO+THT)  *  (RO+REL)  )  )  ; 

%Estimate  initial  angle  of  elevation  using  4/3  earth 

beta43=. 75*betat;  %i.e.  beta (effearth) = (re/re (eff) ) *betat 

rc43=l . 25*rc;  %rc  modified  for  eff  earth 

r043=rc43+SEL;  %R0  modified  for  eff  earthd 

trg43=sqrt  (  {r043+REL)  ''2+  (r043+THT)  •^2-2*  (r043+REL)  *  .  .  . 

(r043+THT) *cos (beta43) ) ; 
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if  trg43==0  %if  trg43  happens  to  be  zero  it  will  cause  a  ?/0  error 
%unless  we  modify  it  slightly 
trg43=le-12 ; 

end  %end  of  check  for  trg43=0 

aO=acos ( (r043+THT) *sin (beta43) /trg43) ;  %Note:  aO  is  always  positive 
if  (r043+THT) < ( (r043+REL) /cos (beta43) ) , aO=-aO; , end  %determine  whether 
%  takeoff  angle  is  +  or  -  using  Pythagorean  theorem  (see  notes, step  5) 

if  -isreal (aO) , a0=le-12; end  %if  aO  is  close  to  zero,  it  may  be 
calculated 

%as  imaginary  or  complex,  hence  this  fix 

%Determine  whether  target  and  radar  are  in  a  duct  together 
if  DUCTFLG>0 

if  DUCTFLG==2 

Mgrad= (pM (2 ) -pM ( 1 ) ) / (ph (2 ) -ph (1)); 
botductht=ph  (2)  +  (pM(3)  •“pM(2)  )  /Mgrad; 
else 

botductht=0 ; 
end  %end  of  DUCTFLG=-2 

if  REL<ph (3) &REL>botductht&THT<ph (3) &THT>botductht 
bothinduct=l ; 
else 

bothinduct=0 ; 

end  %end  of  height  comparison 
else 

bothinduct=0 ; 

end  %end  of  target-radar  in  duct  together  check 

^  ^  ic  ^  ^  -k  -k  -k  -k  -k  -k  -k  ^  ^  ^  ^  -k  -k  -k  -k  ^  -k  -k  -k 

%Initialize  counters, flags, etc . 

hittarg=0; 

hithole=0 ; 

a0toolo=~99 ; 

a0toohi=99; 

a0magtoolo=0; 

a01astgood=99; 

ERRMSG= [ ] ; 
hitbetae= [ ] ; 
stoptrace=0; 
a01oB==99  ; 
a0hiB=99; 
tracecount=0; 
cant find-0; 

while  '-hittarg&-'hithole&-cantf ind  %BEGINNING  of  estimate/iterate  loop 
tracecount=tracecount+l; 

%disp ([’ trace  #  ’  num2str ( tracecount )  ’  a0=’  num2str (aO, 16) ] ) 

al=aO;  %initialize  el  angle  #1 
dir=sign (aO ) ; 

lvl=REL/dh+l-dir;  %initialize  level  to  one  before  REL  (taking  into 
%  account  the  increment  at  the  start ) --Note :  level  1  is 

%  at  the  surface  (there  is  no  level  0) 

betae=0 ; 

BETAES=[0];  %Initialize  beta  accumulator 

TRGAPP-0; 

dBda0=0; 

cosia2=0 ; 

heights= [REL] ;  %Initialize  height  accumulator 
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maxminf lg=0 / 
hitdirt f lg=0 ; 
stoptrace=0 ; 
nuinmms=0  ; 
numhits=0 ; 
hitbetae= [ ] ; 
hitdBdaO= [ ] ; 
stepstohit=  [  ] 

g. 

o 


%Initialize  flag  for  stopping  trace 
%Initialize  max/min  counter 
%Initialize  target-height-reached  counter 
%Initialize  target-height-reached  subtense  storage 
%Initialize  target-height-reached  dBdaO  storage 
%Initialize  target-height-reached  #-steps-it-took- 
to-get-there  storage 


outoduct=0 ; 

while  ~stoptrace  %BEGINNING  of  trace  loop 
lvl=lvl+dir ; 

if  lvl==2  &  dir==-l  %if  level  is  currently  1  above  ground. . . 
hitdirtf lg=l ; 


end 

cosial=cosia2 ;  %for  use  if  there’s  an  extrema 
if  maxminf lg==l  %if  the  previous  point  was  an  extrema... 

%!!!Note:  Determine  direction  by  sign  of  M-gradient — FIX — p4B- 


dir=-dir;  %swap  directions 

Nl=Nm;  %current  N  and  ht  are  at  the  extrema 
hl=hm; 

N2=N (Ivl+dir) ;  Isince  dir’s  been  swapped  around,  this’ll 
h2=h (Ivl+dir ) ;  %  bring  us  back  to  the  last  level  we  were  at 


cosia2=cosia0;  %el  angle  of  ray  is  opposite  that  of  the 
%  level  hit  before  the  extrema 

maxminf lg=-l;  %set  flag  indicating  current  level  is  max/min 
else 

Nl=N(lvl) ; 
hl=h(lvl) ; 

N2=N (Ivl+dir) ; 
h2=h ( Ivl+dir) ; 

%disp([’hl=’  num2str(hl)  ’  h2==’  num2str(h2)  ’  cos(al)=’ 
num2str (cos (al)  ,  16)  ]  ) 

cosia2-(l+ {N1-N2) * (le-6) - (dir*dh/ (RO+hl) ) ) *cos (al) ; 

%disp ( [ ’ cosia2= ’  num2str {cosia2,16)  ]  ) 
end  %end  of  assign-N-and-h-values  routine 

if  cosia2>l  %max  or  min  —  find  ht  of  extrema,  see  Abel, pp . 23 , 24 
A-(l/ (RO+hl) ) * (N2-N1) *le-6/ (h2-hl) ; 

%  A=(l/ra) * (dn/dz) 

B=.5*(l/(R0+hl)+(N2-Nl)*le-6/(h2-hl) ) ; 

%  B=(l/2) (1/a+dn/dz) 

C-l-cosial;  %C=l-cos (alfa) 

delz=dir*abs ( (abs (B) -sqrt (B^2-A*C) ) /A) ;  %ht  difference 
hm=hl+delz; 

Nm=Nl- (hl-hm) * (N1-N2) / (hl-h2) ; 

N2=Nm; 

h2=hm; 

cosia2^1;  %since  el  ang  at  extrema  is  0,  cos(ang)  is  1 
cosiaO=cosial ; 

maxminflg=l;  %set  flag  indicating  next  level  is  max/min 
nummms=nummms+l ;  %count  number  of  max/mins 
end  %end  of  what-if-there ' s-an-extrema?  routine 
if  hl==f loor (hi ) , disp ( [num2str (hi )  ’km  reached’ ]), end 
nl=l+Nl*le-6; 
n2=l+N2*le-6; 

a2=dir*acos (cosia2 ) ;  %is  SIGN  necessary 

psil=2* (nl-n2) / (tan (al) +tan {a2 ) ) ;  %calc  bending  in  layer 
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betal=abs (psil+a2-al ) ;  %calc  subtense  across  this  layer 
betae==betae+betal;  %accumulate  total  subtense  to  this  point 

%disp([^al=*  num2str (al, 16)  '  a2=’  num2str(a2)  ’  psil=' 

num2str (psil, 16)  *  betae=’  num2str (betae, 16) ] ) 

BETAES=cat (2 , BETAES, betae) ;  %store  all  subtenses  for  later 
heights=cat { 2, heights , h2 ) ;  %store  all  heights  for  later 

%Calc  geometric  range  over  layer  and  accumulate 

trgappl=sqrt  (  (RO+hl)  ^2+  (R0+h2)  '^2-2*(R0+hl)*  (R0+h2)  *cos  (betal)  )  ; 

TRGAPP=TRGAPP+t rgappl ; 

%Calculate  how  beta  changed  with  aO  this  layer  (see  Abel, p. 29) 
if  maxminflg==l  %next  level  is  max/min 

gam= (n2-nl ) / (h2-hl ) ;  %intermediate  var  --  n-gradient  in  layer 
dBldaO=-tan  (aO)  /tan  (al )  +  (2*n0*R0’^sin  (aO )  )  /  .  ,  . 

(tan  (al)  *  (n2/gam+R0+h2 )  )“(psil*tan(aO)  )/(sin(al)  )""2; 
elseif  maxminf lg==-l  %current  level  is  max/min 

gam= {n2-nl) / (h2-hl) ;  %intermediate  var  —  n-gradient  in  layer 
maxminflg=0;  %reset  maxminflg 

%Note:  dBldaO  is  the  same  --  symmetric  about  the  extrema 
else  %neither  level  is  max/min 

dBldaO=-tan (aO) /tan (al ) ttan (aO) /tan (a2 ) - . . . 

(psil^tan (aO) / (tan (al) ttan (a2) ) ) * . . . 

(1/ (cos (al) ^sin (al) ) +1/ (cos (a2) *sin (a2) ) ) ;  %Calc  beta 
%  gradient  for  this  layer 

end  %end  of  calculate-dBldaO  routine 

dBdaO=dBdaO+dBldaO ;  %accumulate  total  beta  change  with  aO 

%record  data  every  time  target  height  is  reached 
if  abs (h2-THT)<le-10 

hitbetae=cat (2, hitbetae, betae) ;  Isubtenses  of  "hit"  points 
hitdBdaO=cat (2 , hitdBdaO, dBdaO) ;  %dif ferential  of  "hit"  points 
stepstohit=cat (2, stepstohit, length (BETAES) ) ;  %number  of  steps 
%it  took  to  get  to  each  "hit"  point 

end 

if  DUCTFLG 
if 

REL<ph (3 ) &REL>botductht& (h2>ph (3) | h2<botductht ) , outoduct=l ; , end  %Check 
to  see  if  trace  started  in  duct  else 

%Check  to  see  if  the  trace  has  left  the  duct 

end 

whether  raytrace  needs  to  be  stopped****** 
if  -bothinduct  %radar  and  target  not  in  duct  together 
if  a0>=0  %initial  el  angle  positive  (or  zero) 
if  numinms=-l  I  hitdirtf  Ig  %unwanted  ducting 
stoptrace=l ; 

if  aO>aOtoolo, aOtoolo=aO; ,  end 
elseif  h2==THT  %target  height  reached 
stoptrace=2 ; 
a01astgood=a0 ; 

end  %end  of  el-ang-pos  options 
else  %initial  el  angle  negative 

if  nummms>l  %unwanted  ducting 
stoptrace=3; 

if  a0<a0toohi&DUCTFLG=2 ,  aOtoohi=aO ; ,  end  %if  we’re 

dealing  with  an 

%elevated  duct,  establish  a  new  upper  limit  for  the  trace 

angle 
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elseif  THT>=REL&length (hitbetae) ==1  %target  higher  than 

radar 

%  &  target  height  reached 

stoptrace=4 ; 
a01astgood=a0  ; 

elseif  THT<REL  Itarget  lower  than  radar 

if  (length (hitbetae) ==2) | (length (hitbetae) ==1& . . . 
hitdirt f Ig&abs (hitbetae-betat ) < . 5*hitbetae) 
%target  height  reached  twice  OR  height  reached  once 
%near  the  target  and  trace  hits  the  ground 
stoptrace=5 ; 
a01astgood=a0; 

elseif  ^length  (hitbetae)  &nuinmms==l  %inax/min  reached^  but 
%target  height  never  reached 
stoptrace=8 ; 

if  aO<aOtoohi, aOtoohi=aO ; end 

end 

end  %end  of  el-ang-neg  options 
end  %end  of  el-ang-pos-or-neg  routine 
else  %radar  and  target  ARE  in  duct  together 

if  hitbetae  %If  target  height  has  been  reached 
if  '-outoduct&hitbetae  (length  (hitbetae)  )  >betat 
%radar  and  target  in  duct  together,  ray  goes  past  target 
stoptrace=6; 
a01astgood=a0 ; 
elseif  outoduct 

if  a0>=0  %tgt  ht  reached, trace  left  duct,aO  positive 
stoptrace=2 ; 
else  %a0<0 

if  THT<REL&  ( length  (hitbetae ) —2  |  .  .  . 

(length  (hitbetae)  ===l&hitdirtflg&  .  .  . 
abs (hitbetae-betat ) < . 5* hitbetae) ) 

%target  height  reached  twice  OR  height  reached  once 
%near  the  target  and  trace  hits  the  ground 
stoptrace=5 ; 
a01astgood=a0 ; 
elseif  THT>REL 
stoptrace=4 ; 
a01astgood=a0; 

end 

end 

elseif  a0>0&hitdirtf Ig  %trace  reaches  target  height  on 
upward  swing, then 

%reaches  max,  swings  down  and  hits  the  ground 

stoptrace=6; 

a01astgood=a0 ; 

end 

elseif  nuitmims>l  I  (numiTiins==l&hit  dirt  fig) 

%If  target  height  has  not  been  reached  but  two  extremas 

%have  been  hit,  or  one  &  trace  hit  dirt 

stoptrace=9; 

end 

end  %end  of  radar&targ-both-in-duct-or-not  check 
if  hitdirtf lg&a0<0&'-stoptrace  %initial  el  ang  neg,ray  hits 
%surface,  and  no  other  stop  situation  exists 
stoptrace=7 ; 

if  aO>aOtoolo, aOtoolo=aO ;  ,  end 

end 
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al=a2;  %Prepare  to  jump  to  next  level 
end  %END  of  trace  loop 


%!!!!!!!!!!!!  I  I  ! tEMP !!!!!!!!!!!!!!!!!!! 

%plot  bent  path 

%plotit=input ( ’ Plot  trace?  {y=l/n=0) ’ ) ; 
plotit=0;  %set  plotit=l  to  plot  every  trace 
if  plotit 

gndrngs-BETAES^RO* , 53996;  %conversion  factor  changes  km  to  Nmi 
htsft=heights*3280 . 83;  %conv  factor  changes  km  to  feet 
plot (gndrngs ^ htsft , ’y’ ) 

text (gndrngs (length (gndrngs) ) , htsft (length (htsft ) )  ,  . . . 
num2str ( tracecount ) ,  *  HorizontalAlignment  * ,  *  center  * ,  . . . 

’ VerticalAlignment ^ , ’bottom') 
hold  on 

plot (betat*R0*.53996,THT^3280.83,  '*b'  ) 
disp ( [ ' hitbetae='  num2str (hitbetae)  ]  ) 
disp ( [ ' betat= '  num2str (betat ) ] ) 
disp ( [ ' stoptrace^'  num2str (stoptrace) ] ) 

%disp ( [ 'hitdBdaO='  num2str (hitdBdaO) ] ) 

disp ('Hit  any  key...') 

pause 

end 

%!!!!!!!!!!!!!! END  TEMP  TEST !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 


fQj^  target  reached/inaccessible’^*'^*'^'^^"^^^^"^^"^^ 
%Check  for  target-reached  or  target-inaccessible  based  on  reason 
%for  stopping  trace  (stoptrace) 

%stoptrace=l , 3, 7 , 8 , 9 :  no  check; reestimate  directly 
if  stoptrace==2  |  stoptrace==4  |  stoptrace=-5  |  stoptrace—6 

[dummy, clostsub] =min (abs (hitbetae-betat) ) ;  %determine  which 
%  "hit"  subtense  is  nearest  the  target 

missby=hitbetae (clostsub) -betat ;  %determine  proximity  of  "hit" 
if  sign (missby) ==-l  %store  takeoff  angles  of  nearest  hits 
a01oB=a0 ; 
else 

aOhiB=aO ; 

end 

if  abs (missby) <=angtol, hittarg^clostsub; , end  %determine  if  trace 
%came  near  enough  to  target  to  stop  tracing 

end 

%Flag  to  stop  if  estimation/iteration  isn't  getting  anywhere 
if  '-hittarg&tracecount>maxnumtraces  %If  more  than  20  traces  have 
occurred 

%and  still  the  target  is  not  reached... 

cantfind=l;  %set  flag  indicating  estimator  can't  find  target 

end 

^★★^★^★★★★★★★★★Estimate  new  initial  takeoff  angle********^'*'*^'^'^*'*' 
%Estimate  new  initial  elevation  angle  based  on  reason  for 
%stopping  the  trace  (stoptrace) 
if  -hit targ& stoptrace 
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a01ast=a0;  %Keep  track  of  the  last  aO  (for  various  reasons) 
switch  stoptrace 

case  1  %pos  el  ang,  unwanted  ducting 

if  a01astgood===99  %If  every  trace  has  been  in  duct... 

a0=1.5*a0;  %double  initial  el  angle  to  break  out  of  duct 
else  %If  a  previous  trace  was  out  of  the  duct... 
aO=aOtoolo+ ( aOlastgood-aOtoolo) /2;  %split 
%difference  between  inti  angs  to  brk  out  w/out  jumping 
%over  last  good  trace 

if  abs (a0~a01ast) <=angtol,  hithole=l, end 

end 

case  2  %pos  el  ang,  target  height  reached 

aO=aO+ (betat-betae ) /dBdaO ;  %use  standard  estimator 
case  3  %neg  el  ang,  unwanted  ducting 

if  a01astgood===99  %If  every  trace  has  been  in  duct... 
if  DUCTFLG==1  %i.e.  surface  duct 

a0=-1.5*a0;  %increase  initial  angle  and  make  it 

positive 

%to  break  out  of  the  duct 
else  %i.e.  DUCTFLG=2,  elevated  duct 

a0=1.5*a0;  lincrease  initial  el  angle  to  break  out  of 

duct 


end 

else  %If  a  previous  trace  was  out  of  the  duct... 
aO=aOtoohi+ (aOlastgood-aOtoohi ) /2; 

%split  difference  between  inti  angs  to  brk  out  w/out 

%jumping  over  last  good  trace 

if  abs (aO-aOlast ) <=angtol,  hithole=l , end 

end 

case  4  %neg  el  ang,targ  ht  (higher  than  radar)  reached 
aO~aO+ (betat-betae) /dBdaO;  %use  standard  estimator 
%if  betae>betat&a0<a01ast, a0=. 75*a01ast , end  %if  ducting 
%causes  estimator  to  fail  and  go  the  wrong  direction, 

%fudge  it  down.  Eventually,  betae  will  drop  below  betat 
%and  aOloB  and  aOhiB  will  kick  in  and  guide  it  in 
case  5  %neg  el  ang,targ  ht  (lower  than  radar)  reached  twice 
aO=aO+ (betat-hitbetae (clostsub) ) /hitdBdaO (clostsub) ; 
case  6  %ducting, trace  reached  targ  ht  at  point  beyond  targ  or 

trace 

%reached  target  height  and  hit  the  dirt 
aO=aO+ (betat-hitbetae (clostsub) ) /hitdBdaO (clostsub) ; 
if  abs (aO) <=aOmagtoolo, a 0=1 . 2*a0magtoolo; , end 
%  use  standard  estimator  with  values  at  hit  closest  to 
%  the  target 

case  7  %neg  el  ang,  hit  ground  (in  duct  or  not) 

if  a01astgood==99  %If  every  trace  has  hit  the  dirt... 
if  THT>=REL 


a0=-.25*a0;  %force  pos  el  angle  to  converge  on  soln 
else  %i.e.  Target  is  lower  than  radar 

a0=.75*a0;  %raise  angle  a  bit  in  order  to  eventually 
%reach  the  target  height  while  remaining  negative 

end 

else  %If  a  previous  trace  stayed  away  from  dirt... 
a0=a0toolo+ ( aOlastgood-aOtoolo) /2; 

%split  difference  between  lowest  airborne  ray  and 

%highest  "underground"  ray 

if  abs (aO-aOlast ) <=angtol,  hithole=2, end 

end 

case  8  %neg  el  ang,  tgt  below  rdr,  but  didn't  reach  low  enough 
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then 


if  a01astgooci==99 

a0=1.2*a0;  %increase  angle  to  get  down  to  tgt  height, 

%  regular  estimation  will  do  the  trick 

else 

aO=aOtoohi+ ( aOlastgood-aOtoohi ) /2; 

%split  difference  between  inti  angs  to  brk  out  w/out 
% jumping  over  last  good  trace 
if  abs (aO-aOlast ) <=angtol,  hithole=l, end 

end 

case  9  %rdr&tgt  in  duct  together,  tgt  height  not  reached 
aOmagtoolo=abs (aO) ; 

a0”1.2'^a0;  %arbitary  increment  amount  —  no  sign  change 

end 

if  stoptrace 

%Ensure  estimate  doesn’t  go  beyond  hi/lo  bounds  established  to 
%avoid  hitting  the  ground  and  getting  trapped  in  a  duct 
if  aO<=aOtoolo 

aO=aOtoolo+ (aOlast-aOtoolo) /2 ;  %split  diff  w/last  aO 

estimate 

end 

if  aO>=aOtoohi 

aO=aOtoohi+ ( aOlast^aOtoohi ) /2 ;  %split  diff  w/last  aO 

estimate 

end 

%Ensure  estimate  doesn’t  go  beyond  hi/lo  bounds  established  to 
%prevent  waffling  around  the  target  without  hitting  it 
if  a01oB'-=99&a0hiB~=99Sc  (abs  (a0~a01oB)  >abs  (aOhiB-aOloB)  |  .  .  . 

abs (aO-aOhiB) >abs (aOhiB-aOloB) ) 

%If  the  limits  are  established  &  the  new  estimate  is  outside 

them 

a0= (aOloB+aOhiB) /2 ;  %redefine  aO  halfway  in-between  limits 

end 

%Ensure  estimate  does  not  go  positive  if  target  is  below  radar 

and 

%they  are  not  in  the  duct  together 
if  aO>0&THT<REL&~bothinduct , a0=. 25*a01ast ; , end 

if  hitbetae>100’^betat  %If  initial  takeoff  angle  is  too 
%small,  a  huge  betae  will  result.  This  knocks  the  estimate 
%out  of  that  range  (eventually) 
aO^lO^aOlast ; 

end 


end 


end 

end  %END  of  estimate  and  iterate  loop 

%*-^***^***^^^*** ★Calculate  final  information************** 
if  hithole  %If  target  is  in  a  radio  hole  or  over  the  radio  horizon. . . 
switch  hithole 

case  1  %radio  hole 

ERRMSG= [’ Target  appears  to  be  in  a  radio  hole’]; 
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case  2  %over  radio  horizon 

ERRMSG= [* Target  appears  to  be  over  radio  horizon']; 

end 

elseif  cantf ind&~bothinduct 

ERRMSG= [' Unable  to  trace  to  target,  probably  due  to  ambiguities 
caused  by  ducting.']; 
elseif  cantf ind&bothinduct 

ERRMSG=[ 'Unable  to  trace  to  target.  Radar  and  Target  are  in  the  duct 
together  resulting  in  ambiguous  estimation']; 
else  %Target  reached 
ALFO=aO; 

HEIGHTS=heights (1 : stepstohit (hittarg) ) ; %Only  keep  those  steps 
BETAES=BETAES (1 ; stepstohit (hittarg) ) ; %necessary  to  get  to  target 
%Calculate  geometric  elevation  angle  (i.e.  Inital  angle  for  straight 
%line  path  from  radar  to  target, 
if  TRG==0  %prevents  potential  ?/0  error 
TRG=le-12; 

end 

alfga=asin  {  (RO+THT)  *sin  (betat) /TRG) -"pi/2;  %calc  geo  el  angle  for  long 

and 

alfgb=“alfga;  %short  RO ' s  (see  notes) 

rta=sqrt ( (RO+REL) ^2+TRG^2-2* (RO+REL) *TRG*cos (alfga+pi/2)  )  ;  %calc  two 
rt ' s  to 

rtb=sqrt  (  (RO+REL)  ^2+TRG"^2-2*  (RO+REL)  ^TRG^cos  (alfgb+pi/2)  )  ;  %compare 
if  abs (rta” (RO+THT) ) <abs (rtb- (RO+THT) )  %whichever  rta/b  equals  rt 
indicates 

ALFG=alfga;  %which  alfga/b  is  the  correct  geometric 

else  %elevation  angle 

ALFG=alfgb; 

end 

%Calculate  measured  (apparent)  target  height 

THTAPP=sqrt  (  (RO+REL)  ^2+TRGAPP'^2-2*  (RO+REL)  ^TRGAPP^cos  (ALFO+pi/2)  )-R0; 

%Calculate  height  error,  range  error,  %  error 
THTERR=THTAPP-THT;  %height  error 
THTERR100=THTERR*100/THT;  %height  percent  error 
TRGERR=TRGAPP-TRG;  %range  error 
TRGERR100=TRGERR*100/TRG;  %range  percent  error 

%Save  duct  parameters  for  plotting  later 
if  DUCTFLG>0,DUCTHTS=[botductht,ph(2) ,ph(3) ] ; , end 


end 


A-39 


■  Appendix  B 

Marsden  Square  Numbering  System  for  the  World 


Bibliography 


84th  Radar  Evaluation  Squadron.  Class  Handout,  Radar  Evaluation  Technical  Training, 
Radar  Theory  Vol  1.  84  RADES,  Hill  AFB  UT  July  1998. 

Abel,  Michael  D.  et.  al.  “The  Theory  and  Use  of  a  Raytracing  Model  Developed  at 
USAFETAC.”  USAFETAC/TN-82/005.  USAF  Environmental  Technical 
Applications  Center,  Scott  AFB  Bu,  September  1982. 

Bean,  B.R.  And  E.J.  Dutton.  Radio  Meteorology.  National  Bureau  of  Standards 
Monograph  92.  Washington:  GPO,  1  March  1966. 

Bean,  B.R.  and  G.D.  Thayer.  “Central  Radio  Propagation  Laboratory  Exponential 
Reference  Atmosphere,  “  Jour.  Res.  NBS,  Vol.  63D,  No.  3:  315-317  (June  1959). 

Blake,  Lamont  C.  Radar  Range  Performance  Analysis.  Norwood  MA:  Artech  House, 
Inc.,  1986. 

California  Department  of  Water  Resources,  Division  of  Local  Assistance.  “CIMIS  ETo 
Calculations.”  Equations  excerpted  from  an  article  by  Pruitt  and  Doorenbos  in  the 
Proceedings  of  International  Round  Table  Conference  on  “Evapotranspiration,” 
Budapest,  Hungary,  1977,  n.  pag.  http://wwwdla.water.ca.gov/cimis/cimis/ 
hq/etocal.htm.  3  Aug  1998. 

Hitney,  H.V.  and  J.H.  Richter.  Integrated  Refractive  Effects  Prediction  System  (IREPS), 
Nav.  Eng.  Jour.,  88,  2:  257-262  (April  1976) 

Kalnay,  E.,  et.  al.  “The  NCEP/NCAR  40- Year  Reanalysis  Project,”  Bulletin  of  the 
American  Meteorological  Society,  Vol.  77,  No.  3:  437-471,  (March  1996). 

Kerr,  Donald  E.  Propagation  of  Short  Radio  Waves.  New  York:  McGraw-Hill  Book 
Co.,  Inc.,  1951. 

Ko,  H.W.  et.  al.  “Anomalous  Propagation  and  Radar  Coverage  Through  Inhomogeneous 
Atmospheres,”  AGARD,  CP346:25-1  -  25-11  (1994). 

Ko,  H.  W.,  Sari,  J.  W.,  and  Skura,  J.  P.,  “Anomalous  Microwave  Propagation  Through 
Atmospheric  Ducts,”  Johns  Hopkins  APL  Tech.  Dig.  4,  12-26  (1983). 

Livingston,  Donald  C.  The  Physics  of  Microwave  Propagation.  Englewood  Cliffs  NJ: 
Prentice-Hall,  Inc.,  1970. 


BIB-1 


Patterson,  W.L.  Advanced  Refractive  Effects  Predictive  System  (EREPS)  Version  1.0 
User’s  Manual.  Technical  Document  3028.  San  Diego:  Space  and  Naval  Warfare 
Systems  Center,  April  1998. 

Patterson,  W.L.  Ducting  Climatology  Summary  (DCS)  software  (Version  2.0,  21 

January  1992)  read  file.  San  Diego:  Space  and  Naval  Warfare  Systems  Center,  1993. 

Patterson,  W.L.,  et.  al.  Engineer’s  Refractive  Effects  Predictive  System  (EREPS) 

Version  3.0.  Technical  Document  2648.  San  Diego:  Naval  Command,  Control  and 
Ocean  Surveillance  Center,  RDT&E  Division,  May  1994. 

Patterson,  W.  L.  Historical  Electromagnetic  Propagation  Condition  Database 

Description.  Technical  Document  1 149.  San  Diego:  Naval  Ocean  Systems  Center, 
September  1987. 

Ryan,  Frank  J.  Analysis  of  Electromagnetic  Propagation  Over  Variable  Terrain  Using 
the  Parabolic  Wave  Equation.  TR1453.  San  Diego:  Naval  Ocean  Systems  Center, 
1991. 

Schelleng,  J.C.,  C.R.  Burrows,  and  E.B.  Ferrell.  “Ultra-Short-Wave  Propagation.”  Proc. 
/RE  21,  no.3:  427-63  (March  1933). 

Science  Applications  International  Corporation.  Investigation  of  Refraction  Issues. 

TR  DO  30-01-04.  SAIC,  30  March  1996. 

Smith,  E.K.  and  S.  Weintraub.  “The  Constants  in  the  Equation  for  Atmospheric 
Refractive  Index  at  Radio  Frequencies,”  Proc.  IRE  41 :  1035-1037  (Aug  1953). 

Squires,  Michael  F.  “Height-Error  Analysis  for  the  FA  A- Air  Force  Replacement  Radar 
Program  (FARR).”  USAFETAC/PR-91/014.  USAF  Environmental  Technical 
Applications  Center,  Scott  AFB  IL,  August  1991. 

Weisbrod,  S.  and  L.J.  Anderson.  “Simple  Methods  for  Computing  Tropospheric  and 
Ionospheric  Refractive  Effects  on  Radio  Waves,”  Proc.  IRE  ##:  ??-??  (Sep  1958). 

(Note:  this  entry  has  DTIC  #  ADA  056800) 


BIB-2 


Vita 


Capt  Todd  Pittman  was  born  18  September  1967  at  Scott  Air  Force  Base,  Illinois. 
In  1985  he  graduated  from  Chittenango  High  School  in  Chittenango,  New  York  and 
entered  Grove  City  College,  Grove  City,  Pennsylvania  with  a  four-year  ROTC 
scholarship  to  study  Electrical  Engineering.  Upon  graduating  with  a  Bachelor  of  Science 
in  Engineering  (Electrical  Engineering)  Degree,  Capt  Pittman  was  commissioned  into  the 
United  States  Air  Force. 

Captain  Pittman’s  first  assignment  was  at  the  Air  Force  Communications 
Command  Headquarters  at  Scott  AFB  where  he  worked  as  an  air  traffic  control  radar 
evaluation  engineer  and  a  communications  test  engineer  until  March  1994.  He  then 
traveled  to  Wright-Patterson  AFB  and  Keesler  AFB  where  he  was  assigned  to  the  Air 
Force’s  Communications  Systems  Evaluation  School  as  the  Advanced  Narrowband  Radio 
Course  Director.  While  at  Keesler  AFB,  he  was  selected  for  a  regular  commission.  In 
July  1997,  he  entered  the  Graduate  Electrical  Engineering  Program,  School  of 
Engineering,  Air  Force  Institute  of  Technology.  Upon  graduation,  he  will  be  assigned  to 
the  84*''  Radar  Evaluation  Squadron  at  Hill  AFB,  Utah. 

Permanent  Address:  5888  Fieldstone  Drive 
Cazenovia,  NY  13035 


V-1 


